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We present a detailed analytical study of ultra-relativistic neutrinos in cosmological perturbation 
theory and of the observable signatures of inhomogeneities in the cosmic neutrino background. We 
note that a modification of perturbation variables that removes all the time derivatives of scalar 
gravitational potentials from the dynamical equations simplifies their solution notably. The used 
perturbations of particle number per coordinate, not proper, volume are generally constant on 
superhorizon scales. In real space an analytical analysis can be extended beyond fluids to neutrinos. 

The faster cosmological expansion due to the neutrino background changes the acoustic and 
damping angular scales of the cosmic microwave background (CMB). But we find that equivalent 
changes can be produced by varying other standard parameters, including the primordial helium 
abundance. The low-Z integrated Sachs- Wolfe effect is also not sensitive to neutrinos. However, 
the gravity of neutrino perturbations suppresses the CMB acoustic peaks for the multipoles with 
/ > 200 while it enhances the amplitude of matter fluctuations on these scales. In addition, the 
perturbations of relativistic neutrinos generate a unique phase shift of the CMB acoustic oscillations 
that for adiabatic initial conditions cannot be caused by any other standard physics. The origin 
of the shift is traced to neutrino free-streaming velocity exceeding the sound speed of the photon- 
baryon plasma. We find that from a high resolution, low noise instrument such as CMBPOL the 
effective number of light neutrino species can be determined with an accuracy of a{N v ) ~ 0.05 to 
0.09, depending on the constraints on the helium abundance. 



I. INTRODUCTION 

Neutrinos play a significant role in the evolution of 
the early universe. They arc expected to provide around 
40%, e.g. Ref. [1], of the total energy density during 
the radiation era. The gravitational potentials (metric 
perturbations) induced by inhomogeneities in the pho- 
ton and neutrino backgrounds are comparable. Due to 
the internal photon and neutrino dynamics the poten- 
tials decay when the growing acoustic, for photons, or 
particle, for neutrinos, horizon of the universe becomes 
of the order of the perturbation scale, i.e. as the pertur- 
bation modes "enter the horizon" . This decay, in con- 
trast to the constancy of the potential generated during 
matter domination by freely collapsing cold dark mat- 
ter (CDM), leads to substantial difference between the 
amplitude of the acoustic oscillations in the cosmic mi- 
crowave background (CMB) on the scales that enter the 
horizon before and after the matter-radiation equality. 
For example, in the model without neutrinos the ampli- 
tude generated by equal primordial power on the smaller 
scales is 5 times larger, Ref. [2]. In addition, the grav- 
ity of both the photon and neutrino perturbations at the 
horizon entry boosts CDM peculiar velocities, contribut- 
ing to matter clustering. 

Neutrino contribution to the radiation energy den- 
sity reduces the redshift of the transition from radiation 
to matter domination, bringing the transition closer to 
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CMB decoupling. This too leads to important conse- 
quences to both CMB anisotropies and matter clustering. 
The reasons are the larger amplitude of the acoustic os- 
cillations entering the horizon in the radiation universe, 
larger early-ISW effect from the transition proximity, and 
suppressed growth of matter fluctuations in the radiation 
epoch. But these and other discussed later effects caused 
by neutrino background speeding up the cosmological ex- 
pansion can generally be mimicked by variations of other 
standard cosmological parameters. For example, the red- 
shift of the radiation-matter equality could be reduced 
not by neutrino background but by CDM density being 
smaller than derived from the fits assuming the standard 
neutrino content. 

The internal dynamics of neutrino perturbations bears 
almost no resemblance to the more familiar acous- 
tic physics of the photon-baryon fluid. All the three 
main distinctions below arise from neutrinos being fully 
decoupled and free-streaming since a very early red- 
shift z^dec ~ 10 10 , long before the hydrogen recombi- 
nation and CMB photon decoupling at z 7 dcc — 1090. 

First, the tightly coupled photon-electron-baryon fluid 
supports compressional acoustic waves. These waves are 
little attenuated until the recombination. Neutrino per- 
turbations propagate differently, by means of free stream- 
ing. Neutrinos escape overdense regions in every direc- 
tion; the projection of their velocity on the density gra- 
dient spans the entire interval [—1, 1] (in units c = 1.) 
The dispersion of the perturbation transfer velocity along 
the density gradient, called "directional dispersion" [3], 
damps subhorizon neutrino perturbations inversely pro- 
portional to time. This damping was noted three decades 
ago [4]. But it was quickly realized [5] that, regardless 
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of their evolution, the subhorizon neutrino perturbations 
exert negligible gravitational effects on other species. 

Second, neutrino stress is locally anisotropic. Accord- 
ing to the Einstein's equations, the stress sources the per- 
turbations of space-time metric. The anisotropic stress 
leads to richer structure of the metric perturbations than 
locally isotropic fluids can provide. 

Third, neutrino perturbations propagate with the 
speed of light, exceeding the sound speed of acoustic per- 
turbations in the photon-baryon fluid. As a result, the 
gravitational effect of neutrino perturbations on CMB, 
viewed in real space, extends beyond the acoustic horizon 
of primordial inhomogencities. We find that this leads to 
a unique phase shift of the CMB mode oscillations in the 
presence of neutrino gravity. 

What new physics can be revealed by the imprint of 
neutrino gravity on the more easily observable species, 
such as CMB photons or non-rclativistic matter? The 
considered neutrino signatures probe the ratio of neu- 
trino and photon energy densities, evaluated when the 
observed scales enter the horizon. Complimentary con- 
straints on the universe composition in the radiation era 
are set by the predictions of Big Bang nucleosynthe- 
sis (BBN) for the primordial abundances of light ele- 
ments. The baryon to photon ratios inferred from BBN 
and CMB are in good agreement with each other, but 
the presently low observational estimates of the primor- 
dial 4 He abundance [6-10] favor the effective number of 
neutrinos, N v , at BBN below its standard value 3.04. 
Joint analyses [11, 12] of the current data on the pri- 
mordial 4 He and D/H abundances and of the cosmolog- 
ical constraints considered by WMAP team [13] (CMB 
+ LSS + Lyman a, fit by ACDM model) give the 2a 
limits 1.7 < N v < 3.0. If a neutrino chemical potential, 
characterizing vjv asymmetries 1 , is treated as a free pa- 
rameter to be marginalized over, the limits relax [18] to 
— 1.7 < N u < 4.1. However, the constraints from BBN 
and CMB should be combined with caution. 

One of the reasons is that the probed redshifts are 
separated by many orders of magnitude. The processes 
that determine the BBN yield of light elements extend 
from the freeze-out of z/^, v Tl and (shortly after) v e in- 
teractions at Zydcc ~ 10 10 to the fusion of light nuclei 
at z ns ~4x 10 8 . On the other hand, the CMB multi- 
poles up to I ~ 3000 probe the neutrino density in the 
redshift range from entry ~ 6xl0 4 to z cq ~ 3.2 xlO 3 (as- 
suming the "standard" cosmological parameters [13, 19].) 
Either the photon entropy or the number of uncoupled 
relativistic species per comoving volume may change 2 in 
the considerable span of the universe history from BBN 



1 Any initial differences among the individual u/u asymmetries for 
the 3 generations of active neutrinos are equilibrated by neutrino 
oscillations by the time BBN begins, Refs. [14-17]. 

2 The change of the photon entropy density is tightly constrained 
below 2 7ch cmoq ~ (2 - 5) X 10 6 , Refs. [20-22], by the Planckian 
shape of the CMB spectrum from COBE [23] . Although energy 



to the redshifts probed by CMB. The responsible physical 
mechanisms could be, though are not limited to, heating 
from decays of massive particles or fields, e.g. [27-30], 
such as expected in thermal inflation [31, 32], or cooling 
by thermal contact with other species, Refs. [28, 33]. 

Another reason is that both BBN and CMB constraints 
depend on certain properties of the uncoupled relativis- 
tic species, in addition to their total energy density. 
For BBN, the relevant characteristics include the asym- 
metry between the active neutrinos and their antipar- 
ticles, their interaction and mixing with other species 
beyond the standard model, and the cosmological ex- 
pansion rate, which may be affected by the density 
of other uncoupled relativistic particles - "dark radia- 
tion" (right-handed neutrinos, Goldstone bosons, mod- 
uli, etc.) and by more exotic phenomena such as early 
quintessence, non-minimally coupled fields, or extra- 
dimensional physics. On the other hand, the CMB 
anisotropics and matter clustering do not discriminate 
between active neutrinos, their antiparticles, and other 
relativistic degrees of freedom. But the dynamics of cos- 
mological perturbations in the unseen relativistic back- 
ground becomes important. In this paper we focus on the 
signatures caused by ultra-rclativistic decoupled parti- 
cles. Their energy spectrum need not be thermal. Given 
initial conditions, their gravitational impact on the "vis- 
ible" species indeed depends only on their combined en- 
ergy density, parameterized by the effective number of 
neutrino species N v , eq. (45). 

While the impact of neutrinos on the light element 
production at BBN has been studied in detail, the neu- 
trino features in the CMB spectra are less well estab- 
lished. Their comprehensive analysis and the investi- 
gation of their potential for probing the primordial ra- 
diation of non-electromagnetic origin are presented in 
this paper. The motivations for exploring these fea- 
tures and the related constraints independently from the 
BBN physics include: verification of the "standard BBN" 
model (SBBN); guidance in resolving the tensions be- 
tween SBBN predictions and observational estimates of 
the light element abundances (the tension presently ex- 
ists for He but in future is also conceivable for other 
elements as potentially more sensitive experiments with 
less studied systematics appear); probing the param- 
eter space of extended BBN models in the directions 



released into the photon gas at smaller redshifts can still be redis- 
tributed among the photons by Compton scattering, the photon 
production rates, from double Compton scattering (e-y — > 677) 
and bremsstrahlung (eN — > eN'y), become insufficient to change 
the photon number to its new equilibrium value. This would 
lead to a Bose-Einstein CMB frequency spectrum with a non- 
zero chemical potential; see reviews in Refs. [24, 25] . The present 
agreement between the BBN- and CMB-derived ratio of baryons 
to photons is an additional evidence against a large change of the 
comoving photon entropy density. Of course, the above consid- 
erations do not limit the change of the energy of the uncoupled 
relativistic species; see e.g. Ref. [26] for specific scenarios. 
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of degeneracies {e.g. the degeneracy in the v chemi- 
cal potential - N v plane [18]); constraining the models 
of high energy physics, frequently leading to decoupled 
relics, non-standard BBN, particle decays during or after 
BBN, or modified cosmological expansion; finally, clar- 
ifying robustness of the constraints derived from CMB 
anisotropics and matter clustering. 

The possibility of identifying the background of de- 
coupled ultra- relativistic species with CMB, sometimes 
complemented by other cosmological probes, has been 
analyzed extensively [12, 34-39] in the past using numer- 
ical calculations with Boltzmann integrator codes, such 
as CMBFAST or CAMB/CosmoMC [40-42], or with sim- 
pler codes in that Boltzmann hierarchy is truncated at 
the quadrupole order in a way that mimics free steam- 
ing [43, 44]. Some of this work forecasted future con- 
straints on the density of relativistic species using likeli- 
hood (Fisher matrix) analysis of specific experiments. 

The authors of Refs. [2, 45] noted that the CMB 
modes entering the horizon in the radiation era are 
perturbed less and the CDM modes more in a model 
with a larger neutrino to photon ratio. Later work [44] 
stressed the essential role of neutrino perturbations in 
breaking the degeneracies between N v and the density 
of non-relativistic matter, set by uj m = il m h 2 , cither 
of which affects the redshift of radiation-matter equality 
z cq + 1 = Pm,ol '(p~i,o + Pv,o)- Degeneracies between the 
variation of neutrino density and of other cosmological 
parameters were studied numerically in Ref. [34]. This 
work pointed out that, with a fixed z cq and fixed angular 
size of the acoustic CMB horizon, the remaining CMB 
spectrum variation with N v to the third acoustic peak 
can be practically removed by a same sign change in the 
scalar spectral index n s , and that the matter power spec- 
trum breaks this degeneracy. However, because of the 
normalization of the calculational results by the height 
of the first acoustic peak, the neutrino-induced suppres- 
sion of CMB anisotropy on small scales was explained as 
increased ISW contribution on large scales. This inter- 
pretation propagated into several later papers. We argue 
at the end of Sec. Ill and Sec. V A that this interpretation 
is incorrect. 

The Fisher matrix likelihood analysis of Ref. [34] 
showed that, prior to WMAP, N v could not be con- 
strained by CMB alone. With the WMAP data [46] 
new analyses [12, 35, 38] set the upper limits N y < 7 
or somewhat better if matter clustering or HST data are 
included. Ref. [38] reported a lower limit N v > 1.6 with 
95% confidence level from WMAP only and N v > 1.9 
with HST data added. We find that these constraints 
can be improved dramatically with the future experi- 
ments and become comparable and tighter than those 
presently derived using the standard BBN model from 
the primordial element abundances. 

Recently, Ref. [47] considered the interaction of neu- 
trino perturbations with tensor gravitational waves. The 
problem was reduced to an integro-differential equation 
using the so-called line-of-sight solution for free streaming 



particles, derived previously in a context of photons [48- 
50]. Numerical integration of this equation showed that 
neutrinos suppress the amplitude of the gravitational 
waves entering the horizon in the radiation era and of 
the related S-mode of CMB polarization by about 20%. 
Even on the largest angular scales the neutrino damping 
of the tensor correlation functions is predicted to be close 
to 10%. 

In this work we focus on the more significant and, as of 
now, the only accessible to observations scalar 3 pertur- 
bations. We use an analytic approach. It provides the 
physical insight into the cosmological role of neutrinos 
and helps find a quantitatively small but unique signa- 
ture of neutrino perturbations, the phase shift, which 
turns out to play the primary role in measuring the neu- 
trino background density with CMB experiments. The 
analytical methods developed in this paper are easily ap- 
plicable to the tensor sector and give results consistent 
with Ref. [47]. 

A real-space view of cosmological perturbation dynam- 
ics will be essential for obtaining analytic results for neu- 
trino perturbations, which can not be modeled by a fluid. 
Many equations governing perturbation dynamics in the 
radiation era are integrated trivially in real space. This 
permits analytic calculations that would seem hopeless 
in momentum space. A real space analysis of cosmologi- 
cal perturbations was attempted earlier in Ref. [51] and 
applied to CMB anisotropy in Refs. [52, 53]. We follow 
the plane wave formalism developed in Refs. [54-56]. 

The rest of the paper is organized as follows. In Sec. II 
we introduce a slight modification to the classical defi- 
nition of cosmological perturbation variables. A conse- 
quence is substantial simplification of the evolution equa- 
tions, both for their later solution and conceptual under- 
standing. In Sec. Ill we set up the notations and evolu- 
tion equations for the radiation-matter universe around 
the time of CMB decoupling. Then we study the im- 
pact of neutrinos on the evolution of superhorizon per- 
turbations. In Sec. IV we review the Green's function 
formalism and apply it to find how neutrino perturba- 
tions influence the CMB and CDM modes that enter the 
horizon in the radiation era. A reader not interested in 
the specifics of the analytic calculations can look at their 
results in the figures of Sees. Ill and IV and proceed to- 
ward the discussion of Sec. V. In Sec. V we analyse the 
neutrino signatures in the CMB and matter perturbation 
spectra and either their robustness or degeneracy to the 
variation of other cosmological parameters. In Sec. VI we 
estimate the accuracy of constraining the effective num- 
ber of neutrino species from some planned or proposed 
CMB experiments. We conclude in Sec. VII. 



3 As customary in cosmology, the term "scalar perturbation" de- 
notes the invariance of the perturbation Fourier modes with re- 
spect to the little rotational group: the axial rotations that do 
not change a mode wave vector k. 
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Appendix A reviews the linear cosmological perturba- 
tion theory and summarizes the properties of the used 
metric gauges. In Appendix B we prove that all the mat- 
ter or metric Green's functions in the Newtonian or syn- 
chronous gauges vanish for growing adiabatic perturba- 
tions beyond the particle horizon. Appendix C contains 
technical calculations for Sec. IV. 

All the following formulas imply the metric signature 
(—1,1,1,1). Greek indices range from to 3; Latin from 
1 to 3. Dots denote the derivatives with respect to con- 
formal time dr = dt/a, where a is the cosmological scale 
factor. The universe expansion rate with respect to con- 
formal time is denoted by Ti. = a/ a = aH, where H(z) is 
the proper Hubble expansion rate. 

II. DYNAMICAL PERTURBATION VARIABLES 

In this paper we use predominantly the conformal 
Newtonian, later "Newtonian" , gauge [57, 58] and pa- 
rameterize scalar metric perturbations as 

ds 2 = a 2 (r) [(-1 - 2$)dr 2 + (1 - 2*)dr 2 ] . (1) 

The potential $ determines the gravitational acceleration 
of free-falling objects g = — V$. characterizes the per- 
turbations of spatial curvature in this gauge. (This choice 
of potentials agrees with Ref. [58] . It is related to other 
frequently cited works as: <&thcrc = ^hcrc, * there = $here 
for Refs. [57, 59], and * t here = $here, $thcrc = -fhorc for 
Refs. [2, 45, 60].) In the presence of anisotropic stress, 
provided by neutrino perturbations, the potentials $ and 
^ differ from each other. 

Occasionally, we invoke the synchronous, spatially flat, 
comoving, and uniform density gauges. Their definitions 
and the relations between various metric gauges are sum- 
marized in part 3 of Appendix A. 

A. Coordinate particle number densities 

It appears very useful to describe perturbations of mat- 
ter 4 species in terms of the variables that rate of change 
does not depend on the time derivatives of other per- 
turbation variables. This is not so, for example, for 
the usually considered proper energy density enhance- 
ment 5pa/pa in the Newtonian gauge, nor the proper 
phase-space density perturbation 5f a , nor the local CMB 
temperature anisotropy 6T(h)/T, nor the "effective tem- 
perature" perturbation ST/T + <!>, for all of which the 
corresponding conservation equations involve W or $ + ^ 
terms. Via the Poisson equation, eq. (A46) in Ap- 
pendix A, these terms, which are dominant on horizon 



4 In Sec. II "matter" refers to all the dynamical degrees of freedom, 
whether in relativistic or non-relativistic particles or fields, in 
contrast to non-dynamical scalar metric perturbations. 



scales, bring the time derivatives of other matter pertur- 
bations. This complicates the description of perturbation 
evolution. 

Instead, we characterize density perturbation of 
species a in the Newtonian gauge by a variable 

rf a E^ = ^-3*, (2) 

^a, coo 

where 

S a = faa - prop = (3) 

^-a, prop Pa i Pa 

is the proper particle number overdensity 5 . The latter is 
related to the energy momentum tensor T£ v by eq. (A19) 
in Appendix A. Unless noted otherwise, we suppose that 
the matter species in any group a do not interact non- 
gravitationally with the species of the other groups; hence 
Tg v is well defined and covariantly conserved. Examples 
of the species groups a arc: photon-baryon plasma, neu- 
trinos, or cold dark matter (CDM). 

The species mean velocity and anisotropic stress, de- 
fined with eqs. (A19), for scalar perturbations can be de- 
scribed by a velocity potential u a as Vi a = — Viii a and by 
an anisotropic stress potential 7r a , eqs. (A19-A21). From 
the energy and momentum conservation equations (A41) 
and d a definition (2), 

d a = V 2 u a , (4) 

U a = C 2 a d a - XaU a + V 2 7T a + $ + 3 C 2 * , (5) 

where 6 c 2 = dp a /dp a is the "sound speed" and 

Xa = H (1 - 3c 2 ) (6) 

is the Hubble drag rate for the species a. Eqs. (4) and (5) 
can be combined into a single second order equation 

d a + Xad a - c 2 a V 2 d a - V 4 tt = V 2 ($ + 3c 2 *) . (7) 

For the special cases of CDM and tightly coupled photon 
fluid with negligible baryon density eq. (7) reduces to 

d c + Hd c = V 2 $ (CDM) ; (8) 

d\ - \ V 2 d 7 = V 2 ($ + *) (photon fluid) . (9) 

o 

The variables d a can be interpreted in the Newto- 
nian gauge as the perturbations of the conserved par- 
ticle number densities with respect to the coordinate dif- 
ferential volume, e? 3 r, rather than the proper volume, 



5 If n is the density of any conserved number, its change in a locally 
inertial frame for a closed volume V equals dn/n = — dV/V = 
dp/(p + p), given the energy conservation d(pV) + pdV = 0. 

6 Except for the generalized proof of superhorizon conservations 
in Sec. II B, wc will assume that for all the species a the local 
pressure p a is uniquely specified by the local energy density p a - 
This assumption is general enough to apply to photon-baryon 
plasma, cold matter, massless or massive neutrinos, and constant 
vacuum energy. It is not valid for a classical field (quintessence) 
or modified Hubble expansion (Cardassian energy.) Without this 
assumption eq. (4) should be replaced by eq. (39). 
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a 3 (l — 3^)d 3 r. The change of these densities with time 
is determined only by the particle flux into a unit coor- 
dinate volume, cq. (4). ft does not explicitly depend on 
the metric evolution, affecting the proper number den- 
sities S a - More formally, d a corresponds to the gauge 
invariant quantity 



d a = S a - 'ZH L 



V 2 x 



(10) 



where Hl and x parameterize the general perturbation 
of spatial metric gij as given by eqs. (A25-A26). There- 
fore, the variables d a coincide with the particle number 
overdensities 5 a in the gauge where the spatial metric 
is unperturbed - the "spatially flat" gauge of Sec. A3. 
The density perturbations d a are simply related to the 
"conserved curvature" perturbations of Ref. [61] on the 
hypcrsurfaces of uniform energy density of species a, 

(uniform, a) 



as d a — 3£ a . 



Likewise, we eliminate \& from the equations of per- 
turbation dynamics in phase space. For this purpose we 
define a variable df a (r, r, q, n): 



df 

df a =6f a + q^V 
aq 



(11) 



where 8f a (T,T,q,n) is the perturbation of the proper 
phase space distribution for species a, f a {r,q) is their 
unperturbed density in phase space, and the particle co- 
moving momentum q = qh is defined in Appendix A by 
eq. (A53). The linearized Boltzmann equation (A57) in 
terms of the variable (11) reads 

(df a )'+ - ntfi (df a ) = q ^ mVi ( £ ¥ + 1 *V (12) 
e dq \t q J 

Summation over the repeated index i = 1, 2, 3 is implied. 
Having in mind the applications to collisionless particles, 
we dropped the collision term and the terms involving 
the time derivatives of f a and of df a /dq. 

The number of phase-space coordinates in the Boltz- 
mann equation can be reduced by one, Ref. [62], when 
the mass of the particles is negligible relatively to their 
average kinetic energy, so that e ~ p and the f a perturba- 
tions propagate with the same speed, the speed of light, 
regardless of the particle energy. Defining a function 



D a (T,r, n 



_ 3 J q 2 dqqdf a {Y,q,n,T) 



(13) 



4 Jq 2 dqqf a (q) 
and integrating both sides of eq. (12) over q 3 dq, we find 



D a + mViDa = -3niVi(* + *) 



(14) 



The variable D a (r, r, n) is related to the energy-averaged 
phase space distribution perturbation of Refs. [48-50] as 
F a {r,T, n) = (4/3)L> a + 4*. If the free-streaming par- 
ticles had the thermal velocity distribution at their de- 
coupling then the temperature perturbation of the parti- 
cles moving in a specified direction n is 5T a (r, r, n)/T a = 
D a /3 + *. 



Eq. (14) is formally solved by 

D a (r, r,n) = D a , m (r - nr, n) - (15) 
-Zn^J\r' (* + *)| r ,, r _ A(r _ T , ) . 

In the following subsection we will show that for adia- 
batic initial conditions Z? a .i n (r, n) is independent of n. 
It will be related to the conserved superhorizon value of 
the spatial curvature perturbation ( in the uniform den- 
sity gauge (the "Bardeen's curvature" [63]) as 



D a . ;„ (r, n) = -3Cin(r) 



(16) 



for all a. 

Any scalar perturbation D a (r, r, n) can be described 
by scalar multipolc potentials {d; )a (r, r)};=o.i,... as 



D a (n) = ^(-l) i (2/ + l)P i 



V l di, a , (17) 



where Pi are Legendre polynomials. Since Pi(fJ-) contains 
only the powers y}~ 2q with q = 0, 1, . . . \ l/2\ , the gradi- 
ent operator enters the right hand side of eq. (17) only 
through natural powers of n^Vi or V 2 . The potentials 
di a are gauge invariant for I > 2. 
As follows from eqs. (A55, 11, 13), 



(18) 



where n = 1, hq = —1, n % = n^, and ()„ stands for 
/ Gi 2 f2n/47r. Substituting the multipole expansion (17), 
remembering the definition of the variables d a , u a , 
and 7r a , eqs. (2, A19-A21), and using that for the ultra- 
relativistic particles p a = p a /3, we find 



d ,a = d a , di. a = U a , d 2 ,a = ^ ■ (19) 



d>i, a dynamics follows from eqs. (14, 17) and the relation 
fi p, j i_ 

— n+i + 21+1 



di.a ^ + 1 dl-l,„ , 



V'd l+ha + S n + (20) 



(The Kronecker symbol Sn in the last term should not 
to be confused with a density perturbation.) One can 
write a formal integral solution of these equations by ex- 
panding the line-of-sight solution (15) over the spherical 
harmonics. 7 



7 Namely, 
dl, a (r, r) = 3 



CinW + 



(21) 
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The scalar metric perturbations $ and \& are deter- 
mined from the linearized Einstein equations, eqs. (A42- 
A45). In terms of the introduced dynamical variables. 

V 2 *-3 7 * = -/(d + 3Hu) , (23) 
* - $ = 3 7 tt , (24) 

where 7 = A?:Ga 2 (p + p) and 

d = '^x a d a , u = ^2x a u a , 7r = ^x a 7r a , (25) 

a a a 

where x a = (p a + p a )/(p + p)- Eqs. (23-24) for * and 
$ are non-dynamical constraint equations. These con- 
straints, however, do not limit one's freedom of setting 
the initial matter distribution potentials di^ a . 

The above equations of matter dynamics, e.g. 
eqs. (8, 9, 14) for the model with dark matter, photon 
fluid, and neutrinos, that are complemented by the el- 
liptic gravitational equations (23-24) form a well posed 
closed system. They are the basis for our subsequent 
analytical analysis of the perturbation dynamics. 



B. Superhorizon conservation of the coordinate 
number densities 

To illustrate the conservation laws prior to their gen- 
eral derivation, we start from the model where all the 
scalar dynamical (matter) degrees of freedom can be de- 
scribed by scalar potentials d^ a that satisfy the evolution 
equations of the form (4-5,20). Below, A is a character- 
istic comoving scale of perturbation variation in space. 
For a harmonic mode it can be taken as A = 1/fc. Since 
in the subsequent sections we are interested in the uni- 
verse evolution long after the inflation, we first assume 
that the comoving Hubble scale H _1 (r) grows through- 
out all the considered time. Then it is natural to choose 
the zero of time r at the formal limit of the equations 
7Y _1 — > 0. This choice will be implied whenever we refer 
to a specific r value, in the context of post-inflationary 



expansion. All of the above assumptions will be lifted by 
the end of the section. 

We set the initial conditions at a time T m < A as 
di.a ~ ( r in) Z - Such initial conditions are natural for grow- 
ing modes, where perturbations are finite for t — > 0. If 
the global intrinsic curvature K does not dominate the 
Friedmann expansion and the universe does not inflate 
(w > -1/3) then, by eqs. (A15,A13), H ~ 1/r and 
7 ru 1/t 2 . Then for our initial conditions and r <C A 
we can drop all the V 2 terms in the evolution equations 
up to an error 0(r 2 /A 2 ). We thus find 

d a , (26) 
u a ~ c 2 {d a + ZHu a + 3*) -Hu a + <I> , (27) 
di, a — 2T+\ di-i,a (free stream. I > 2) , (28) 

d + 3Hu + 3V ~ , (29) 

* - $ = 37?r . (30) 

The corresponding solutions scale for r in « t < A as 
di.a ~ t 1 and $ ~ \& ~ 1. Specifically, the coordinate 
particle number density perturbations d a are constant, 
up to 0(r 2 /A 2 ) corrections. Adiabaticity of the initial 
perturbations has not been assumed for this result. 

If the initial conditions are adiabatic 8 , from cq. (A24), 
the density and velocity perturbations of all the species 
in the Newtonian gauge are related in the limit r/A — > 0, 
A = const as 

S a = —3Hu a = S (adiabatic, r/A — > 0) , (31) 

where S is a same function of r and r for all the species a. 
Also, by eq. (2), 

d a = d (adiabatic, r/A — > 0) (32) 

with d = S — 3^ for all a. Then the terms in eq. (27) 
that are proportional to c 2 cancel by eq. (29). Hence, all 



+ T df 4Sr + *(*.')]} 

where ji are spherical Bessel functions. This equation is obtained 
from eqs. (15-16) and (17) by noting that for any analytic func- 
tion /(r) 

/(r-nr) = e -™- v -/(r)= (22) 
= f>lV(2i + 1) P ; (?^L\ SM-iVr) /(r) . 

In eq. (21), the operators ji(kr)/k l and j'^kr) /k l_1 , with k 2 — > 
— V 2 , are well defined as their Taylor expansions involve only 
even powers of fe, hence, only integer powers of the Laplace op- 
erator V 2 . If the perturbations in eq. (21) refer to a single spatial 
harmonic plane wave with a wave- vector k then —V 2 does be- 
come k 2 . 



8 We define a perturbation as "adiabatic" (a curvature perturba- 
tion) if in some space-time coordinates all the proper matter dis- 
tributions or fields and their proper rate of change, all smoothed 
over an arbitrary comoving scale A, appear unperturbed in the 
limit r/A — » 0, A = const. (Without the assumption H ~ 1/r, 
the relevant limit is [64] W _1 (r)/A — * 0, A = const.) In this, 
"superhorizon" , limit the space-time metric in the considered 
coordinates remains perturbed. We prove in Appendix B that 
for any wave-vector k there exists a non-decaying perturbation 
mode with such properties. This mode satisfies the conditions 
of adiabaticity of Ref. [65] and Ref. [59]: As shown in this sub- 
section, under very mild natural assumptions about the internal 
matter dynamics, the corresponding comoving curvature pertur- 
bation Ti is constant beyond the horizon, up to 0(t 2 /A 2 ) correc- 
tions. Since after a coordinate change (A23) the proper energy 
density and pressure of species a appear perturbed as 8p a = p a Sr 
and Sp a = pair in the linear order, all the ratios Sp a /p a and 
Sp a /p a in any other gauge in the above limit are equal. 
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the velocity potentials evolve identically on superhorizon 
scales in the leading and next to the leading orders in 
r/A: 



-Hu a + $ 



(adiabatic, r <C A) . (33) 



The same is true for all d a — d, which remain constant 
by eq. (26). The leading order evolution of the I > 2 mul- 
tipoles depends on the internal dynamics of the species. 
For example, these multipoles vanish identically for per- 
fect fluids but they grow for neutrinos. Nevertheless, the 
given definition of adiabaticity demands that d^ a — *• in 
the r — > limit. This is consistent with our previous ob- 
servation that, for the "growing mode" initial conditions, 
eq. (28) yields di, a ~ r l . 

In addition to the coordinate particle number den- 
sity perturbation d, two other perturbation variables are 
known to be constant for adiabatic perturbations on su- 
perhorizon scales. These are the spatial curvature per- 
turbation in the comoving gauge 



TZ = * + Hu , 



(34) 



cqs. (A36, A38), Refs. [59, 65-67], and the curvature per- 
turbation in the uniform density gauge £, eq. (A39) [63]. 
(The sign of 1Z and £ in this paper coincides with the sign 
of the Newtonian potentials. Most of the references use 
the opposite sign.) Since the scalar metric perturbations 
are not dynamical but are fixed by matter perturbations 
with constraint equations, one may expect a simple re- 
lation between the conserved matter perturbation d and 
the metric perturbations 1Z or £. Indeed, the comparison 
of eq. (2) and eq. (A39) gives immediately that 



-3C 



It also follows from eqs. (A40, A17) that 



n = c + o, 



(35) 



(36) 



Thus, up to 0(r 2 /A 2 ) corrections, for the growing adi- 
abatic perturbations and all the species a, 



d a (r «),r) = — 3Cin(r) (adiabatic) 



(37) 



where Cin(r) is the time-independent superhorizon value 
of the curvature perturbation £. Substituting this result 
in eq. (17) and remembering that for I > 1 di. a ~ t 1 — > 0, 
we also find 



D a (r < A,r,n) = -30„(r) (adiabatic) . (38) 

The conservation of D a (ii) or, for non- interacting parti- 
cles with non-negligible mass, of df a (q,h) in eq. (11) is 
also evident from eq. (14) or (12), in which all the terms 
with gradients can be dropped for superhorizon growing 
perturbations by the same arguments as before. 

Eq. (7) assumed that local pressure of the species is 
uniquely determined by their local energy density. With- 
out this assumption, for all the mutually non-interacting 



groups a of matter species or fields, the energy conserva- 
tion T^.fj, = gives 



PaSpg — PaSpa 
(Pa+Pa) 2 



(39) 



(We emphasize that Tg v and all the variables in eq. (39) 
are defined Appendix A independently of the nature of 
the species self- interaction.) As previously, V 2 u a can be 
dropped for growing superhorizon perturbations. Since 
energy density and pressure perturbations transform un- 
der a gauge transformation (A23) as Sp a = Sp a + p a Sr 
and Sp a = Sp a +p a Sr, the additional term on the right 
hand side of eq. (39) is gauge invariant . If for a considered 
group of species a, which may interact among themselves 
but do not couple to the other species, the energy den- 
sity and pressure become homogeneous as r/A — > in 
some coordinate frame, this term is initially zero in any 
frame. We call such initial conditions for the species a 
"internally adiabatic" . For them, all the right hand side 
of eq. (39) vanishes and d a is constant beyond the hori- 
zon. Generalization of the arguments that follow eq. (41) 
gives that d a starts changing only in the order 0(r 2 /A 2 ). 

If all the species a are perturbed internally adiabat- 
ically, which is automatic for single-component perfect 
fluids, then we showed that all the d a are constant for 
growing modes beyond the horizon. However, the vari- 
ables d a need not be equal for different a's if the overall 
perturbation is not adiabatic. In this case 



1 



(40) 



in general changes outside of the horizon as the species 
enthalpy abundances x a , eq. (A10), vary during the ex- 
pansion. This is essentially the curvaton mechanism of 
Ref. [68], see Ref. [69] for a modern version, converting 
isocurvature into curvature perturbations. 

Most generally, any system with locally interacting 
matter and Einstein gravity possesses a covariantly con- 
served energy- momentum tensor T^ v . Therefore, the 
scalar perturbation variables d = S — 3^ and u are always 
well defined with eqs. (A19-A20). From the covariant 



conservation T° M ., 



0. 



d= VV 



pSp — p5p 
(p + p) 2 



(41) 



In Appendix B we show that if a growing adiabatic per- 
turbation is initially localized in a spatial region then all 
the matter and gravitational Newtonian gauge potentials, 
including u, vanish beyond the particle horizon of this re- 
gion. Then, by the Gauss's theorem, the velocity diver- 
gence term in eq. (41) has zero integral over any volume 
enclosing the particle horizon. For the initial conditions 
that are adiabatic as defined in footnote 8, the gauge in- 
variant quantity (inon-ad = (p 3p — P&p) I (p + p) 2 tends 
to zero in the limit r/A — > 0, A = const. If, motivated by 
either the equivalence principle or analyticity of the sys- 
tem dynamics, we accept that the metric perturbations 
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generate d n on-ad only in 0(r/A 2 ) order (in the frame 
where the matter is initially unperturbed but so in any 
other frame because of the gauge invariance of d non -ad) 
then in the leading and next to the leading orders in r/A 
the variables d and £ = —d/3 are constant. So is constant 
the comoving gauge curvature perturbation 1Z, eq. (36), 
up to 0(r 2 /A 2 ) deviations. 

If a stage of inflation, defined as cosmic expansion with 
positive acceleration d 2 a/dt 2 , hence contracting comov- 
ing Hubble scale 7i _1 (r), is also considered, the condition 
r/A -> should be replaced by W" 1 (r)/A -> 0, Ref. [64]. 
Then this limit and the related definition of adiabaticity 
in footnote 8 in the inflating universe apply to the future 
rather than the past. 

Finally, what about the universe evolution that begins 
with an inflationary stage and proceeds to the canonical 
Big Bang? Ref. [65] defines adiabaticity and consid- 
ers the conservation laws in the limit A — > oo. How- 
ever, the initial conditions for modes with finite wave- 
lengths can be different from those for the infinitely large 
scales. To accommodate such physically viable possibil- 
ities, we define the adiabatic initial conditions for any 
fixed finite spatial scale A = 1/fc on the stage of in- 
creasing 7i _1 (r) by formally evolving the perturbation 
with the post-inflationary equations backward in time 
to Hr 1 — > 0. (Similarly, by evolving forward in time to 
JiT 1 — > during the inflation.) This approach allows one 
to quantify the primordial non-adiabaticity (admixture 
of isocurvature modes), which can be probed by CMB 
or matter spectra, at any k. The lesser the primordial 
non-adiabaticity and the "tidal" 0(H.~ 2 / \ 2 ) dynamical 
deviations are, the better the proved conservation laws 
apply. 



Thus in much of this work we will be interested in 
the perturbation dynamics at the redshift z > 10 3 . Bar- 
ring the possibilities of noticeable early quintessence [71] 
or the Cardassian modification of Fricdmann expan- 
sion [72, 73], the background expansion at that time can 
be described by a flat radiation-matter model. The ra- 
diation energy density is provided by CMB photons and 
neutrinos, which mass becomes dynamically relevant only 
at z < m v /(3kT v fi) ~ 200 m, y /(0.1 eV) and is neglected 
here. (The published WMAP [13] 95% CL limit on the 
neutrino masses is m„ < 0.23 eV.) The massive matter 
consists of cold dark matter, c, and baryons, b. 

We begin from establishing notations convenient for 
the radiation-matter model. The linearized gravitational 
equations, e.g. (23-24) or (A42-A45), involve the re- 
duced enthalpy background density 7 = 4wGa 2 (p + p), 
cqs. (A11-A13). For baryons and CDM, with p a = 
Pa,a/o? and negligible pressure, it equals 



7b(c 



(42) 



a 2a 

where H rc { = 100kms _1 Mpc _1 andwft( c ) = Vl b ^h 2 . The 
present WMAP constraints [13] on uj b and u) m = uj b +uj c . 
assuming the standard neutrino content, are ui b = 0.024± 
0.001 and to m = 0.14 ± 0.02. 

The reduced enthalpy of photons, 7 7 = (16/3)7rGa 2 /? 7 , 
at a given redshift 1/a is fixed by today's CMB temper- 
ature T 7 , = 2.725 ± 0.002 K, Ref. [74], as 



7 7 



3a 2 a 2 



(43) 



where w 7 = Q 7 ft 2 w 2.47x 10~ 5 (T 7 , /2.725 K) 4 . For neu- 
trinos, 



III. RADIATION-MATTER UNIVERSE 

If the primordial perturbations are nearly adiabatic, 9 
the inhomogencitics in the neutrino background may af- 
fect only those CMB and CDM perturbations that en- 
tered the horizon while the radiation fraction of the uni- 
verse energy was non-negligible. As for the perturba- 
tion modes with larger wavelengths, the number density 
perturbations d a = — 3fi n remain frozen and the higher 
angular multipoles <iz>i ia of the photon and matter dis- 
tributions negligible until the horizon entry. By the time 
these modes enter the horizon and the species distribu- 
tions start evolving, the neutrino energy density pertur- 
bations are too small to have a gravitational impact on 
their evolution. 



9 The following argument and its conclusion do not apply to isocur- 
vature initial conditions, considered recently for neutrinos in 
Ref. [70]. 



7, = a,7 7 , a v = fj- = 7 - (^) 3 N V ~ 0.23 N v , (44) 

where the standard Big Bang Nucleosynthesis predicts 
N v ~ 3.04, assuming 3 Standard Model neutrino genera- 
tions and zero neutrino chemical potential. The effective 
number of neutrino species N u ^ 3, first, because neutri- 
nos share some of the energy of e + e _ , annihilating soon 
after the neutrino decoupling peak, Ref. [75]. Second, be- 
cause this energy, most of which heats the photons after 
the annihilation, is somewhat reduced by finite tempera- 
ture QED corrections [76, 77] . The physics of both effects 
is concisely reviewed in Ref. [78] . 

Of course, here we allow JV„ to be a free parameter. 
It characterizes the energy density of all the decoupled 
ultra-relativistic species at the considered moment, im- 
plied being after e + e~ annihilation but before CMB de- 
coupling, 



■Ny — Prcl decoup I g (ll) 



(45) 



This parameter may have a non-standard value due to 
cither an unaccounted change of p v or p 7 , or due to the 
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density of additional weakly interacting ultra-rclativistic 
species {X). The latter would presumably decouple 
at very high temperature when the universe contained 
more relativistic degrees of freedom than during the neu- 
trino decoupling. As the particles such as heavy leptons, 
hadrons, W and Z bosons, Higgs fields, superparticles, 
etc. become non-relativistic and annihilate, the entropy 
shared by the coupled photons, electrons, and neutri- 
nos increases. Since the comoving entropy density of the 
decoupled species X remains unchanged, their contribu- 
tion to the parameter (45) may become substantially be- 
low unity. A light field carrying gx effective degrees of 
freedom, with the fermionic ones multiplied by 7/8, that 
decoupled when the remaining particles in thermal con- 
tact had g(T < x dec) degrees of freedom contributes to the 
ratio (45) as 



AN,, 



7 



fix 



g(T <x dec) 

9yev 



(46) 



were g jel , = 43/4. So, for example, a hypothetical neu- 
tral Majorana fermion or a scalar Goldstone boson that 
decoupled when the remaining relativistic degrees of free- 
dom were composed of all the fields of the minimal super- 
symmetric standard model {g(T < x dec) = 915/4) would 
give AN U ~ 1.7 x 1CT 2 and AN V ~ 9.7 x 10~ 3 corre- 
spondingly. 

It is sometimes convenient to use the ratio of comoving 
time t to a characteristic time of the radiation-matter 
energy equality, and the ratio of scale factor a to its value 
at the equality: 



where 



(1 + a v ) w 7 



U)„ 



(47) 



(48) 



3.5- 10 3 V 1.69 



1 



0.3 x 0.7 2 



and 



Note for reference that r = 4r e (l— r)/r, dr — —4T e dr/r 2 



and 



where 



4a cq (l — r)/r . In terms of : 



n 



my 



Pr 



1 



av 



(52) 



(53) 



is the neutrino fraction of the total radiation energy den- 
sity p r = Pl + Pu - R v m 0.408 for N v = 3.04. The for- 
mulas describing the supcrhorizon perturbation modes, 
Sec. IIIC, become very compact if the mode evolution is 
parameterized by r. For example, see eq. (75) for the well 
known k = growing mode of the gravitational poten- 
tial in the radiation- matter neutrinolcss model, Ref. [60] . 
The radiation and matter domination limits of these for- 
mulas are easily read off by setting r to 1 and corre- 
spondingly. 



A. Perturbations in the radiation era 

When the universe energy density is dominated by pho- 
ton gas and ultra-relativistic neutrinos, for all of which 



w a — cz = tt, then 



i 

T 

(r) 



(54) 



R v , and x^ for 



In the radiation era x 
any non-relativistic species is negligible. The evolution 
equations (9), (14), (23-24) in this regime become 



1 „ 
d 1 — — V dy 

O 

Dy + TliViDv 

t 2 V 2 * -6* 



= V 2 ($ + 1') , (55) 

= -3n 4 V,(* + $) , (56) 

= 2d+ -u , (57) 

T 

= *-^. (58) 



2(V2 - 1) H le{ ' 



130 Mpc 



'l + a v /0.3 x 0.7 2s 



1.69 



(49) 



The Friedmann equation for the radiation-matter uni- 
verse in terms of the variables (47) reads (da/df) 2 = 1+a, 
yielding 



a 



1 2 
T+-T 2 

4 



We find it useful to introduce the variable 
f 1 2 



1 



i + VTT 



a 



(50) 



(51) 



On the scales well inside the acoustic horizon, A <C 
r/y/3, the gravitational terms on the the right hand side 
of eqs. (55, 56) are negligible. The oscillating photon 
acoustic modes and the free-streaming neutrinos decou- 
ple from each other. However, the phase and the ampli- 
tude of the acoustic oscillations is set by the perturbation 
dynamics during the horizon entry, when the gravity of 
neutrino perturbations played a significant role. Like- 
wise, the neutrino perturbations affect the dark matter 
peculiar velocities, evolving according to eq. (8), as the 
matter received a gravitational boost from the radiation 
when the fluctuations entered the horizon in the radia- 
tion era. 

If the neutrino density is negligible (i?„ — > 0) then 
the above equations have compact analytic solutions in 
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either Fourier space, eqs. (110,117,127), or real space, 
eqs. (109,116,126). When the neutrino gravity is ap- 
preciable, it appears difficult to track the evolution of 
the Fourier modes analytically through the horizon entry. 
But in Sec. IV we succeed with an analytical approach 
in real space. 



B. Perturbations in the matter era 

When massive matter, with pressure and anisotropic 
stress being negligible comparatively with its energy 
density, becomes the dominating component then by 
eqs. (52) and (A13) 



7 



(m) 



6 



(59) 



Provided the energy density perturbations are also dom- 
inated by massive matter, giving negligible n, by cq. (24) 
the gravitational potentials $ and *S> are equal. After 
baryons decouple from CMB photons at Zd ~ 1090, we 
can also take Sp « 0. Then by eq. (A44) 



$ + = 



(60) 



The corresponding non-decaying solutions are time in- 
dependent on all scales. The constant $ and ^ modes 
that enter the particle horizon after Zd are easily related 
to the epoch-independent primordial curvature perturba- 
tion Cm- Indeed, by applying eq. (A38) to superhorizon 
scales, where 1Z = Cm by eq. (A40), one finds 



(m,A»7i- 1 ) _ q,(.m,\>H-*) 



5^ 



(61) 



When $ = ^ = 0, we can rewrite eq. (14) for ultra- 
relativistic neutrinos or for photons after their decoupling 

as 



e 



cir 



off 



0, 



where 



ef(r,r > n) = -X> + * + * 



(62) 



(63) 



(If the free-streaming particles were in local ther- 
mal equilibrium at their decoupling then 0° ff (r, r, n) 
is the relative temperature perturbation of the parti- 
cles propagating in the direction n plus the gravita- 
tional redshift correction: 9° ff = ST a (r, r, n)/T a + $.) 
For superhorizon adiabatic perturbations in the mat- 
ter era Qf (r<A,r,n) = -Cm + $ + * = |*(r), 
eqs. (16,72). The corresponding, "Sachs- Wolfe" [79], so- 
lution of eq. (62) at later time is 



6f(r,r,n) = i$(r-nr) 



(64) 



The related multipole potentials dj „, eq. (17), for a sin- 
gle Fourier harmonic < i ) k(r) =Re(^4e lkr ) follow from 



eqs. (63, 64) and eq. (22) in footnote 7 on p. 5 as 



d a = 



sin(A:T) 



- 6 



U>l,a 



k l 



$1 



(65) 



The evolution equations for the perturbations of other 
species also have simple analytic solutions when the grav- 
itational potentials are time independent. We do not 
write these solutions here because the linear perturbation 
dynamics in the matter era has been thoroughly stud- 
ied in the past, and neutrinos, while relativistic, do not 
modify it. Of course, the power spectra of perturbations 
in the matter era are affected by neutrinos through the 
change of the effective initial conditions for the modes 
that entered the horizon prior to matter domination. We 
discuss this modification of CMB and matter power spec- 
tra in Sec. V. 



C. Superhorizon Scales 

The I > 2 multipoles d^ a of the phase space distribu- 
tions for the non-relativistic CDM and baryons are neg- 
ligible in the linear regime. They are also small for CMB 
photons, isotropized by scattering prior to hydrogen re- 
combination. The integral solution for the multipoles 
of free-streaming particles, eq. (21) in footnote 7, shows 
that for neutrinos on superhorizon scales d/.„ ~ r Cin- 
Particularly, the neutrino anisotropic stress potential 
7r„ = |d 2 , v is °f the order of r 2 Ci n - Then, by eqs. (24, 52), 
^ — $ = 37 I ,7r„ ~ R v r 2 Cin- Hence, in the radiation era, 
when r — ► 1, the neutrino anisotropic stress leads to split- 
ting the Newtonian gauge potentials $ and \? even on 
superhorizon scales, Ref. [57]. The effect disappears in 
the matter and later eras, after the energy-momentum 
tensor of neutrinos becomes negligible comparatively to 
that of non-relativistic species. 

In Sec. II B we showed that for superhorizon growing 
adiabatic perturbations d a = — 3Ci n = const for all the 
matter and radiation species. The superhorizon evolu- 
tion of the I > 1 multipoles di^ a and gravity is described 
by eqs. (33,28-30) of Sec. II B. In that section we also 
observed that for adiabatic perturbations all the velocity 
potentials u a are equal, up to 0(r 2 /A 2 ) corrections, to 
the momentum-averaged velocity potential u. Then, by 
eqs. (28,19), 



717, = U 

15 



(r « A) 



(66) 



Combining this with eqs. (33, 29-30,35) we find a closed 
equation 

4 4 

7T„ + 2HK + -7„7T„ = — Cm (t < A) . (67) 

Its non-decaying numerical solution for N v = 3 is plot- 
ted in Fig. 1 a) with the solid line. Given the neutrino 
anisotropic stress potential 717,, the superhorizon value of 
the velocity potential u can be calculated from eq. (66). 
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FIG. 1: The evolution of superhorizon adiabatic perturbations in the radiation-matter universe, a) Neutrino anisotropic stress 
potential, b) The Newtonian gauge gravitational potentials. On both plots, the solid curves show the full result for 3 neutrino 
species. The dotted curves correspond to neutrino species. The dashed curves on plot b) are the sums of the leading and 
subleading terms in the expansion of the potentials in R„, for 3 neutrino species. The dashed vertical lines show a/a cq at CMB 
decoupling, given the cosmological parameters of Ref. [46] . 



The gravitational potentials $ and \& are then obtained 
from eqs. (33,24). The potentials corresponding to the 
N v = 3 numerical solution of cq. (67) are plotted in 
Fig. 1 b) with the solid lines. 

We observe that, first, coincidentally, the Newtonian 
potential $ is almost unchanged during the radiation- 
matter transition if N„ ks 3. Second, when radiation 
is dynamically significant, the sum $ + 4" is smaller in 
a universe with a larger effective number of neutrinos. 
This sum governs the propagation of CMB photons, as 
seen from eq. (7) with c 2 = i. The following analytical 
analysis quantifies these observations. 

In the radiation era = 1/r, 7,^ = 2R 1/ /t 2 , and 
the growing mode solution of eq. (67) is 



(r,-r<A) 



15 + 4R L 



Hence, from eqs. (66,33,24), 



(r,r«A) _ 



1 



2Ch 
3 



20. 



1 + t^R v 



^(r,r«A) = [\ + ± R v ) $<r,T<X) 



(68) 

(69) 
(70) 



Relation (70) between the potentials in the radiation era 
was previously derived in Ref. [57]. 
In the matter era 

ftM = 2/t and the 7j,7Tj, term in 
eq. (67) is negligible. Then 



_(m, t«A) 



r 2 20, 



25 3 

and wc obtain the conventional result: 



$ (m,r<SA) = ^(m,T«A) _ £ ^ 



(71) 
(72) 



c.f. eq. (61). 

In the intermediate regime cq. (67) has no simple ex- 
act solution. But the physics of the superhorizon pertur- 
bation dynamics in the presence of neutrino anisotropic 
stress can we studied analytically by expanding the solu- 
tion in the powers of R v . The calculations in the zeroth 
and the first orders are straightforward and are given be- 
low. 

In the zeroth order, i.e. when the neutrino fraction R v 
is negligible, the gravitational potentials $ and ^> are 
equal. Then using eqs. (33, 29) and remembering that on 
the superhorizon scales d = — 30n we have 



l(«v— >)' = <„ 



(73) 



This relation is easily integrated when r of eq. (51) is 
taken for the evolution variable: 



30 U 



(74) 



The gravitational potentials then follow from eq. (29) as 

$M = „ M = ^ (9 + ^ + ^ ^ . (75) 

This fluid limit solution, known in more lengthy forms 
before from Ref. [60], is plotted in Fig. lb) with the 
dotted line. The anisotropic stress potential 7r^ of a trace 
amount of neutrinos can be found by the integration of 
eqs. (66, 74) as 



r 0R, 



+ o) = r 2 U(r) 20r 
15 3 



(76) 
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where 



Mr) 



3-r 2 



2r 2 kii 



The function f^(r) is plotted in Fig. 2 a). Its radiation 



to a/a eq ~ 6.0, and has a small negative value at smaller 
redshifts. 

The photon density perturbation d 1 on superhorizon 
scales remains constant and independent of $ and $ evo- 
lution. Thus on the scales entering the particle horizon 



(r — ► 1) and matter (r — ► 0) era limits are 1 and § cor- at the redshifts z < z cq /6 (for w m = 0.14 and 3 neutri- 



respondingly. In Fig. 1 a) we compare the leading order 

solution r Ky R "'~*^ /t 2 , dotted curve, with the previously 
found numerical solution for N v = 3. 

To determine the 0(R U ) terms in the gravitational po- 
tentials, we rewrite eqs. (66,67) as 



1 



71 ( a2 ") = C" 



6K 



(77) 



where we substituted the result (52) for 7„. The 0(R V ) 
solution for the velocity potential is obtained by using 
the O(l) solution (76) for 7r„ on the right hand side of 
cq. (77). After its integration, 



with 



Mr) 



3r^ 
U 



4rCi, 
45 



f u R„ + O(Rl) , 



(78) 



(2 - r)f„ - 1 



Hence, by eqs. (29-30), 



80. 



1 - r 



$ = $(^o)_^ URv + 0{R 2 ) 
45 

45 



(79) 



nos, these scales exceed the acoustic horizon at recom- 
bination threefold or more) the potential variation and 
the induced by it ISW effect in the CMB temperature 
anisotropy are little affected by neutrino perturbations. 
Even more so, the background and the adiabatic per- 
turbations of relativistic species play no role in the late 
ISW effect, caused by the global potential decay during 
the universe transition from matter to dark energy domi- 
nation. By this time, their energy density is dynamically 
irrelevant. 



IV. STUDYING RADIATION ERA WITH 
GREEN'S FUNCTIONS 

The evolution of cosmological perturbations in the lin- 
ear regime may be studied by superposing pcrturbative 
solutions (Green's functions) that are localized in real 
space. It is convenient to consider the Green's functions 
that vary with only one spatial coordinate [56], say x. 
They are related to the Fourier space perturbation modes 
by one-dimensional Fourier transformation. For example, 
for the curvature perturbation £, 



where 



dk 
2^ 



e lkx C(r,k) . 



fa(r) = (2-r)f u 



Mr) 



All the functions f u , /$, and tend to 1 in the radiation 
era limit r — > 1, and to in the matter era limit r — > 0. 
/$(r) and are plotted with solid curves in Fig. 2 b). 

The dashed lines in Fig. 1 b) show the sums of the lead- 
ing and sublcading terms in the analytic solutions (79) 
with R v set to its standard value 0.408, corresponding 
to N v = 3.04. As seen from the plots, the 0{R„) ap- 
proximations describe the main features of the numerical 
solutions rather well. The about 11% smaller than pre- 
dicted splitting between the potentials iff and <& in the 
radiation era corresponds to 4i?„/15 « 11% smaller ac- 
tual value of the anisotropic stress i: v than it is given by 
the leading order formula (76), c.f. Fig. la). 

By eqs. (9, 14), both after the photon decoupling and 
when the baryon loading is negligible prior to the de- 
coupling, the photon dynamics is affected only by the 
sum $ + . This sum depends on the neutrino abun- 
dance as 



Normalizing the Fourier modes to 
C(r^0,fc) = Cin 



(81) 



(82) 



where C in is a fc-independent constant, we see from 
cq. (81) that 



C(t -> 0,2:) = Cin<5 D (x) 



(83) 



=- ¥ PA(r)-/,(r)] 

The combination 2/$ — is plotted in Fig. 2 b) with the 
dash-dotted line. It vanishes at r ~ 0.55, corresponding 



(Sd(x) denotes the Dirac delta function.) Thus the con- 
sidered Green's function describes the linear evolution of 
a sheet-like curvature perturbation that was created on 
the whole plane x = and is independent of the y and z 
coordinates. 

The initial ratios of the perturbations for different 
species should be specified as well. In our analysis the 
initial conditions are assumed adiabatic, but the Green's 
function method can be generalized to incorporate ad- 
mixture of isocurvature perturbations. The right hand 
side of eq. (82) could also be chosen as A(ik) n where 
A is constant and n is assumed natural. The resulting 
Green's functions would also be initially localized, with 
0(R U ) . (80) linv^o £(t, x) = A5^ (x). These initial conditions are 
even or odd with respect to the parity transformation 
x — > —x if n is even or odd correspondingly If the ini- 
tial conditions for the relative species perturbations have 
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FIG. 2: The functions f n , , and the combination 2/<j. — that appear in the 0(R V ) order of superhorizon perturbation 

evolution, as considered in the main text. The evolution variable r is defined by eq. (51). The radiation density domination 
corresponds to r — » 1, the matter density domination to r — > 0. 



the same parity, this parity will be preserved for all r. 
It is convenient to impose even initial conditions, as im- 
plied by eq. (83), for adiabatic perturbations and odd for 
isocurvature ones. 

We discuss only the "growing" mode Green's functions, 
corresponding to growing Fourier modes in eq. (81). The 
decaying solutions of the evolution equations are irrele- 
vant if the primordial perturbations were generated many 
c-folding before the scales of our interest entered the hori- 
zon. 

The following two observations prove extremely handy 
in calculating the Green's function. First, applying the 
inverse Fourier transformation to eq. (81) and setting k 
to 0, one finds the following simple connection between 
the integral of a Green's function over all the space and 
the superhorizon Fourier modes of the same variable: 



/+oo 
dx((T, x) = C(r,fc^0) , (84) 
-oo 



or an analogous relation for any other perturbation vari- 
able. As it was shown in Sec. II B, for adiabatic perturba- 
tions the right hand side of a sum rule such as eq. (84) is 
time independent for the curvature or density perturba- 
tions (, 1Z, d a , D a , or df a . It vanishes for the growing adi- 
abatic perturbations of all the I > multipole potentials. 
In the radiation era, the k — ► limit for any other pertur- 
bation, i.e. involving gravitational potentials, is trivially 
calculable for any neutrino density from the results of 
Sec. IIIC. 

Second, the adiabatic Green's functions with even ini- 
tial conditions (83) identically vanish beyond the parti- 
cle horizon of the original perturbation, \x\ > t, for all 
the considered perturbation variables in the Newtonian 
gauge, including the potentials <& and 'J. This result is 
proven in Appendix B. This is a non-trivial statement, 



taking into account that the / > 1 multipole and grav- 
itational potentials are not locally measurable physical 
quantities and their dynamics is not necessarily causal. 10 



A. Green's functions for phase space distributions 

Neutrinos and decoupled photons are described by 
their distributions f a (r, r, q, n) in the phase space 
(r, q = qri). A scalar perturbation of f a that does not 
depend on y and z coordinates must also be indepen- 
dent of n y and n z : Sf a = 5f a (r,x,q,fJ,), where /i = n x . 
For ultra-relativistic free streaming particles, the energy- 
averaged distribution D a (r, x,(i), eqs. (13,11), satisfies 
the transport equation 

D a +fiVD a = -3 A ;V(* + $) , (85) 

where V = d/dx, c.f. eq. (14). The corresponding mul- 
tipole potentials di M (r,x), eq. (17), equal 

(-l) l V l dt, a = J^Pt(n)D a (v) . (86) 

To illustrate the application of Green's functions to 
neutrino dynamics, we first consider the free streaming of 
masslcss particles in a time- independent gravitational po- 
tential, as it is the case in the matter era, Sec. IIIB. The 
transport equation satisfied by the effective temperature 
Green's function 9f (r,i,/i) = § D a + $ + eq. (62), 



For gravitationally interacting perfect fluids the linear evolution 
of the Newtonian gravitational potentials turns out to be causal, 
as it can be shown by generalizing the formalism of Ref. [56] to 
arbitrary fluids. This is not true in general. As a simple coun- 
terexample, consider an absolutely inelastic collision at x = 
of two identical sheets of ultra-relativistic particles that are or- 
thogonal to the x axis and move toward each other with opposite 
velocities. At the moment of the collision the system anisotropic 
stress at x = disappears. This changes \& — <I> instantly through- 
out all the space. 
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becomes 



0° 



M V0f = 0. 



(87) 



Given the initial conditions 0° ff (0, x, //) = Qin5r>(x), it is 
solved by 



Qf(r, x, fx) = e in S D (x - fir) 



The multipole potentials di_ a corresponding to D a = 
3 (0|f - $ - *) follow immediately from eq. (86). Re- 
membering the definition of the Legendre polynomials 

i riV-iy 



we find 



dla- 



2 l +m 



(- 



6{r - \x\) - 3J;o($ + , (89) 



where is the Heaviside step function. In the matter 
era $ = * = $i n 5v(x), in = §i n /3, and the Fourier 
transform of the above formula reproduces the modes of 
eq. (65). 



B. Neutrino distributions 

In the radiation era the evolution equations (55-58) 
can be converted from partial into ordinary differential 
equations with respect to a dimcnsionlcss variable 



A 



(90) 



Indeed, the growing modes of such perturbations as g? 7 , 
D v , or ^ during radiation domination have the form 
f(r,k) = /(fcr). The corresponding Green's functions, 
c.f. eq. (81), scale as 



$(t, x) = - $(x) , $(t,x) = -V(x) , 

(91) 

d 7 (r,x) = -d y (x) , D„(T,x,fi) = -D v {x,n) ■ 

T T 

All the first or second partial derivatives of a function 
/(t, x) = f{x)/ T are complete \ derivatives, which we 
denote by primes, times some power of r: 



/ = -(Xf'y/r 2 , V/ = f'/r 2 , 
f=(X 2 f)"/r 3 , Wf=-(xf)"/r 3 

V 2 / = f"/T 3 . 



(92) 



The powers of r can be canceled out of all the terms in 
the evolution equations (55-58). For future references, 



let us note that 

f(r,k) 
We define 



-foe 



d\e 



-ikrx 



fix) ■ 



$ + EE 



$ + * 



* - $ 



(93) 



(94) 



The gravitational potential $_ is sourced directly by neu- 
trino anisotropic stress as described by eq. (58). On the 
other hand, the motion of photons in the radiation era 
and of neutrinos, eqs. (55,56), is affected by <I> + only. 

Applying differentiation rules (92) to the neutrino 
transport formula (85), we obtain an easily integrable 
equation 

[( X -M)^]' = 6/i$ / + . (95) 
Since the Green's functions vanish for |x| > 1, 

(x-M)A> = 6//$+ . (96) 
Eq. (96) does not constrain Z)„ at x = It is satisfied by 



A/(x,m) = Pv{v) Mx - m) + 



6/( 

x-m 



with any function p„(/z). (Even when 4> + and so the right 
hand side of eq. (85) are identically zero, there are non- 
zero D v solutions that describe free streaming neutrinos 
in Minkowski space.) 

The function p v (fi) in eq. (97) must be fixed by the ini- 
tial conditions. Considering rfc — > limit of relation (93) 
and remembering eq. (38), for any < 1 we find 



dx D u (x, M) = D„(rk — >• 0, fj,) = -3Ci, 



(98) 



Substituting solution (97) into the left hand side of the 
above identity, 



For the multipoles 
DiAx) = 



6fi 

dx 

-l X ~ V 



<Mx) 



(99) 



1 dfi 



eqs. (97, 99) give: 



DiAx) = -3 



2Cinfl(x) + 



1 , , $+(x')x^(x) + Mx)xW) 

d X ; 

■ l X - X 



(100) 

f(i-W). 



In Fig. 3 a) the solid line shows the neutrino density per- 
turbation Dq^(x) = d u (x) that is obtained from this 
equation and the potential $ + in the limit R v — > 0, 
eq. (109), when the integrals in eq. (100) are easily taken. 



15 



/(a 



(x - a) r 



-, 71=1,2,. 



x — a 
In \x — a\ 



(x - a)"* 1 

In \x — a 



71=1,2, 



n(x — a) r 

In |x — a 

In |ar — o| + 
n (i — a) n 

- In 2 \x - a\ 



TABLE I: The integrals of singular generalized functions. 
Additional valid formulas are obtained by simultaneously 
multiplying the expressions on the left and on the right by 
sign(a- — a). 



Tabic I gives the definite integrals of several singular 
functions, interpreted as generalized functions. This ta- 
ble is easily understood by noting that the generalized 
functions corresponding to the expressions on the left 
are defined as the derivatives of the less singular expres- 
sions on the right. Of course, the generalized integration 
agrees with the conventional one on any interval on that 
the Ricmann integral exists. 

Table II lists the Fourier transforms of singular gener- 
alized functions. It is useful to remember that the Fourier 
image of an even real function is even and real and of an 
odd real function is odd and imaginary. 



D. Gravitational potentials 



C. A note on generalized functions 



Now we turn to the linearized Einstein equations and 
solve them consistently with the dynamical equations for 
the relativistic matter. First, we differentiate eq. (58) 
twice to obtain 



The expressions under the integrals in eq. (99) or (100) 
are singular at x' = X- The value of the integrals depends 
on how the singularity is treated during the integration. 
Physically, this ambiguity corresponds to resolving the 
last, divergent, term in eq. (97) outside of an interval 
X G [n — £i , P + € 2\ an d approximating the Di „ structure 
inside the interval by the first, <5-function term in eq. (97). 
The most direct approach is to take the integral Cauchy's 
principal values, implying e± = ti — > 0. 

Soon we will encounter the integrals of even more di- 
vergent expressions, such ass -2 or x ~ 2 In x, for that even 
the Cauchy's principal value does not exist. Neverthe- 
less, we can proceed with their meaningful calculation 
if all the singular expressions are understood as gen- 
eralized functions. A detailed mathematical treatment 
of the latter can be found in Ref. [80]. The physical 
meaning of these calculations is clarified by the follow- 
ing two theorems. The first one states that any gen- 
eralized function is a finite order generalized derivative 
of a continuous function. 11 The second that differentia- 
tion of a generalized function f(x) multiplies its Fourier 
components f(k) by ik. Green's functions can be for- 
mally defined as Fourier integrals of perturbation modes, 
eq. (81). Generalized functions provide a consistent, ele- 
gant formalism for their manipulation even when the in- 
tegrals diverge in the Riemann's sense. One could avoid 
the divergences by working only with sufficiently smooth 
potentials of singular real space perturbations, e.g., con- 
sidering $ instead of dp oc V 2 $, on small scales. But 
the use of generalized Green's functions simplifies and 
streamlines the calculations. 



11 To be precise, any finite order derivative of a continuous func- 
tion f(x) with | /| bounded as \x\ — » oo by a finite power of |x| 
defines a generalized function and, conversely, any generalized 
function can be presented as such a derivative. 



= 2R V D 2 , V 



(101) 



where we applied the last of eqs. (19). The perturba- 
tion i?2 jl /(x) on the right hand side is given in terms of 
the potential $ + by eq. (100) with I = 2. 

Second, we note that for adiabatic perturbations in 
the radiation era Sp/5p = 1/3. Then we can eas- 
ily eliminate all the matter perturbations from Einstein 
eqs. (A42, A44-A45) to find 



3 



V 2 * 



V 2 $ 



3* + $ = 



(102) 



Using eqs. (92), we obtain a relation that can be inte- 
grated once trivially, giving 



I)$' + + ( x 2 -l) 



2 X I> + = . (103) 



The general solution of the above equation is 



*+(x) 



const — F- (x) 



d X ' 



x' 2 -i 



(x' 2 



(104) 
(105) 



The integration constant on the right hand side of 
eq. (104) is unambiguously defined for all x if on ly we 
specify how the integral in the second equation is under- 
stood for x > ~ j when the integration path encounters 

singularities at = ±-^=. As discussed in the preceding 
subsection, we treat the singular expression under the 
integral as a generalized function. Then, given a certain 
, one can integrate the singular terms using Table I. 
The constant in eq. (104) may differ among the x inter- 
vals (—oo,—;^~), (^75^75)1 an d (75) °°)- For example, 
if R v = then <l_ = 0, -F-(x) = and l»+(x) vanishes 
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+ 00 dip itp x 



S^ix-a), n = 0,1,2,. 
signx 
(x-a)" n , n=l,2,... 
X - ™ signx, n = 1,2,... 
ln|xl 

X - n ln| X |, rc = l,2,... 
X^ n signx In |x| , n = 1, 2, 



/(<p) = /+^dxe- i ^/(x) 



(*V) n e~*°*' 
2 



(—0221 ,o n_1 

(n-1)! V 3 



sign 95 e 



( (n-i)T y" 1 signy>[ln|y>| -V>("0] 



TABLE II: The Fourier transforms of singular generalized functions. The derivation of these results can be found in Ref. [80]. 
The values of ip(n), the logarithmic derivative of the gamma function, for a natural argument follow recursively from ?/>(n + l) = 
i + ip(n) and ip(l) — —7, where 7 ~ 0.5772 is the Euler's constant. In every case, shifting the transformed function argument 
by a constant a, as displayed in the first and third lines, multiplies the Fourier image by e~ laip . 



for |x| > ^ but not for |x| < c.f. eq. (109). As 
shown in Appendix B, with the initial conditions that 
are adiabatic and satisfy eq. (83), the metric must re- 
main unperturbed beyond the particle horizon |x| = 1- 
From this we immediately conclude that that const = 
in the interval (— 00, — ^=). Taking into account that 

F_(x>l) = because $'_(x) is odd, we see that the 
constant also vanishes for x > 77^- Denoting the value 
of the constant in the interval (— -^=, by p$, we thus 
have 



in tightly coupled radiation fluid: 

dR»->o) ,_ ['siaift cos<#, 



>{T,k)=2& 



(110) 



with ip s = fcr/vo- 

The next iteration step, giving the potentials in 0{R V ) 
order, is performed in Appendix C. Fig. 3 b) shows the 
zcroth order <J> + (dotted) and the 0(R U ) corrected $ + 
(solid) and (dashed) potentials for N v = 3.04. 



P<t>0 



1 

7! 



F-(X) 



. (106) 



p$ may be calculated from the results (69-70) for super- 
horizon Fourier modes in the radiation era that give 



E. Neutrino effect on CMB perturbations 

The Green's function for the photon density pertur- 
bation can be easily found in terms of the gravitational 
potential <J> + . From eqs. (55,92) we obtain 



1 _ 1 

dx$+(x) = - 
1 l 



15^ 



2 i 

3 



=^ . (107) 



In the following subsection we obtain another, equivalent 
but easier to apply, condition fixing p$ . 

The system of the integro-differcntial equations (101) 
and (105-106) for the pair ($+,$'_) can be solved by 
iterations starting from the solution in the limit R v — > 0. 
The latter is 



(fl„-»o) _ 3\/3 Cii 



3~ X 



(108) 
(109) 



as immediately follows from eq. (101) and eqs. (106- 
107). Note that the Fourier transform (93) of the last 
Green's function matches the well known potential modes 



7 



2$ 



+ ■ 



The general solution of this equation is 



ixl " ^= 



2$i 



X 



(111) 



(112) 



The delta function pre-factor p 7 is fixed by the condition 



!_i dxd-yix) = d 7 (kr -> 0) 



Kin, giving 



(113) 



The pre-factor p 7 can be related to the constant p$ , ap- 
pearing in eq. (106), by applying the "Poisson law" (57). 
When eq. (57) is considered as a relation among the 
Green's functions, the only delta-function singularity 
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FIG. 3: a) Adiabatic Green's functions for neutrino (solid) and photon (dashed) number density perturbations in the radiation 
era. The neutrino fraction, R v , of the radiation density is assumed infinitesimal, b) Adiabatic Green's functions for the 
gravitational potentials $± = ± <3?)/2 in the radiation era. The solid and dashed curves are the sums of the 0(R°) and 
0(R V ) terms for three neutrino species. The dotted line is $+ = $ for R v — > 0. 



appearing on its right hand side is the one provided 
by the photon density perturbation (112). As for the 
left hand side, where \f r = $ + + <&_, the only delta- 
function comes from the double derivative of the term 
(X 2 - |) 6 - |x|) in eq. (106). The equality of 
these contributions requires 



p$ = -Vo(l - R v )p 1 



(114) 



Substituting eq. (106) in (113) and eliminating p$ with 
the relation above, we obtain 



1 



Pj = 



\-2R v 



1, _ 
2 W 



d X F-(x) 



(115) 



Calculating p$ from the last two equations is somewhat 
easier than from eq. (107). 

Now we have all the analytic tools to analyze how neu- 
trinos affect CMB perturbations. The evolution of metric 
perturbations without neutrinos is given by eqs. (108- 
109). Then the photon density Green's function follows 
from eqs. (112, 115) as 



4*- 



♦0) 



-30, 



-3&flx 



l 

73 



(116) 



Its Fourier transform (93) leads to the photon density 
Fourier modes in the radiation era: 



d^°\r,k) = -30, 



2 sin ip s 



cos ip 6 



with ip s = kr/\/3. In particular, without neutrinos the 
photon density modes oscillate under the acoustic hori- 
zon (<p s ^> 1) as a pure ip s cosine. 

The predictions for both the phase and the amplitude 
of the photon mode oscillations differ when the gravity 



of neutrino perturbations is taken into account. The os- 
cillations of the Fourier modes on subhorizon scales are 
described by the singular terms in the real space Green's 
functions. For the photon density (112) these are the 
6- function and (x ± -t^) -1 singularities at x = ±- 



l 

V3- 



2r-> 



where 



r 7 = $+(l/V3) 



+ ... , (118) 



(119) 



and the dots stand for more regular terms. The Fourier 
transform of eq. (118) follows from the first and third 
lines of Table II, where n is set to and 1, as 



d~ t (t, k) = 2 (p-y cos tp s 



sm ip s 



O^- 1 ). (120) 



A non-zero phase shift with respect to the cos tp s oscil- 
lations is generated whenever r 1 ^ 0. By eq. (119) this 
can happen for adiabatic perturbations if only some per- 
turbations propagate faster than the sound speed in the 
photon fluid, and thus are able to generate metric pertur- 
bations beyond the acoustic horizon. This is the case for 
the neutrino perturbations, propagating with the speed 
of light, Fig. 3 a). 

The values of p 1 and r 7 in eq. (118) are calculated 
in 0(R V ) order in Appendix C. With its results (C6) 
and (C7), the mode (120) can be presented as 

d 7 (r, k) = 3Ci„(l + A 7 ) cos {ip s + Sip) + O^" 1 ) , (121) 



(117) where 



dip 



-0.2683R v + O(Rl) 
0.1912^ + 0(^2) 



(122) 



As demonstrated in Fig. 4 a), our theoretical predictions 
are in excellent agreement with numerical calculations 
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for the radiation era, at the redshift z = 10 7 , obtained 
with CMBFAST [42]. The numerical calculations show 
that the 0(R 2 ,) corrections contribute to A 7 and Sip less 
than 10% when N v ~ 3. 

Our analytical results (121-122) can be compared with 
numerical fit (B7) of Ref. [2] , giving for the photon over- 
density in the radiation era 

( HS flt) 9^(rfc^ 0)cos^ = SCinCos^ 

2(l + f^) (l + &Rv) 

This formula misses the phase shift due to neutrino per- 
turbations but describes the decrease of the oscillation 
amplitude with the increase of neutrino fraction remark- 
ably well. In the 0(R V ) order it coincides with the re- 
sult (122) by better that 1%. 

Although the primordial magnitude of the cosmologi- 
cal perturbations is unmeasurable directly, one can detect 
the amplitude change of the CMB acoustic oscillations, 
predicted by the first of eqs. (122), by comparing CMB 
and dark matter density fluctuations. The latter, how- 
ever, are themselves affected by neutrinos. In the next 
section we find the leading, 0(RJ), corrections to the 
CDM density perturbation modes entering the horizon 
in the radiation era. 



F. Neutrino effect on CDM perturbations 



(ip s 1), when the potential term in cq. (124) can be 
dropped, the general solution for CDM Fourier modes 
should be of the form 

d c (r, k) = -6Ci„(l + A c ) (imps + 7 - i + Scj + (128) 

The values of the integration constants A c and 8c are de- 
termined by the mode dynamics during the horizon entry 
and are sensitive to the gravity of neutrino perturbations. 
The real space calculations in Appendix C give that 

A c ~ 0.2297i?„ + 0{Rl) , 

(129) 

5c = -0.6323i?„ + 0{Rl) . 

The density perturbation (128-129) in 0(R V ) order is 
compared with the radiation era (z = 10 7 ) CMBFAST 
calculations in Fig. 4b). When kr 3> 1, the 0(R V ) an- 
alytical results underpredict d c variation for N v change 
from to 1 by 11.6% and from to 3 by 23%. Since 
in the second case R v changes by 2.1 times more than 
in the first one, the twofold increase in the relative error 
is consistent with the origin of this deviation from the 
0{R 2 ) corrections. 



The evolution of the CDM coordinate density pertur- 
bation d c is given by eq. (8). In the radiation era it 
becomes 



d c + -d c = V 2 $ . 

T 



(124) 



Using the CDM Green's function ansatz d c (r, x) = 
d c (x)/ T an d the differentiation rules (92), canceling the 
common factors 1 /t 3 , and integrating the resulting equa- 
tion once trivially, we find 



(X 2 dc)' ~ Xdc = 



(125) 



In Appendix C we show that in the R v — > limit and 
with the adiabatic initial conditions this equation gives 
for x + 



In Fourier space 

d c R ^°\r, k) = - 6Ci„ ( In ip a + 7 - \ ci ip s + ^ 
\ 2 p s 



(126) 



(127) 



with ip s = fcr/Vo. 

A finite neutrino fraction of the radiation energy den- 
sity R v affects the gravitational potential on the right 
hand side of eq. (124) and so the matter density pertur- 
bation d c . On the scales well inside the acoustic horizon 



V. NEUTRINO SIGNATURES IN CMB AND 
MATTER SPECTRA 

Decoupled neutrinos affect observable cosmological 
probes both by the gravity of their perturbations and by 
the change of the cosmological expansion rate due to the 
contribution of the neutrino background to the universe 
energy density. The first effect is prominent when cos- 
mological perturbation modes enter the horizon in the 
radiation era. The corresponding modifications of the 
photon and CDM perturbations were found in the previ- 
ous section. The perturbations remain to be propagated 
to the later epochs and related to observable statistical 
power spectra. These tasks are addressed in the current 
section. 



A. CMB power spectra 

Theory overview 

While photons are tightly coupled, their number den- 
sity perturbation d 7 satisfies the equation 



v . KRb ■ 2v7 2 . 
cl-y + — a-y — c s V a 7 

1 + Rb 



(130) 



2r d V 2 d\ + V 2 ( $ 



1 + Rb. 
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FIG. 4: a) Numerically calculated photon number density perturbation d 7 in the radiation era for and 3 neutrino species N v 
(dots) versus the theoretical prediction (117) for N v = (dashed) and its rescaled and phase shifted asymptotic form (121, 122) 
for N u = 3 (solid), b) Similar comparison for the dark matter density perturbation d c and N v = 0,1,3. The theoretical 
predictions are given by eqs. (127-129). In all the cases the 0(R„) terms in the analytical formulas are neglected. 



It follows from eq. (7) applied to the photon-baryon fluid 
with 



Rb(r) 



dp 
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Ref. [45]. The Silk damping term 2r d V 2 d 7 in eq. (130) 
owes its origin to both partial photon diffusion and to the 
lagging of baryons behind photons in imperfectly coupled 
photon-baryon plasma. While these effects arc minor, the 
damping coefficient equals [81] 



T C (T) 



1 - 



14 



1 



15(1 + R b ) (l + R b ) 2 



where 



rc(r) = 



an e a Tho 



(132) 



(133) 



The above result for and eq. (130) are valid while 
t c <C min(l/fc,T). The damping increases substantially, 
Ref. [58], within the CMB last scattering surface, where 
the imperfect fluid approximations fail. 

The general solution of eq. (130) can be obtained for 
subhorizon Fourier modes, fer> 1, using WKB approx- 
imation, Refs. [45, 82]. For the monopole of the photon 



effective temperature perturbation 



u 0,7 — 



(<5T 7 (r, r,n)) a 
T 7 



$ + * , (134) 



(131) it gives 



jyj cos {kS + Sip) - i? 6 $ , (135) 



(l + ifc)V 

where the size of the acoustic horizon S and the Silk 
damping length xs equal 



.At 1 , 4(r) = / T d d T ' . (136) 



S(r) 



The solution (135-136) takes into account that on the 
subhorizon scales photon-baryon fluid and neutrinos con- 
tribute negligibly to the gravitational potentials $ and VP, 
primarily generated by CDM. Hence, these potentials do 
not vary substantially over a single period of acoustic os- 
cillations. It also assumes that Td <C 1/fc, which is a nec- 
essary condition for the validity of eqs. (130) and (132). 

The photon-baryon plasma velocity potential, affecting 
the CMB anisotropy through the Doppler effect, is easily 
found from eq. (4). For Fourier modes, 
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36^ ^ AVi e~ k2x " 
k 2 ~ 



k (l + i? 6 ) 3 / 4 



sin (kS + Sip) , 



(137) 



where the last two equalities arc valid within the approximations that were applied to WKB result (135). 
The present CMB temperature anisotropy observed in the direction n. 



T 



oo I 



Q(il) = J2J2 e i™Yhn(n) , 



(138) 



1=0 m=-l 



can be written as the line of sight integral [42]: 



6(n) = / dr [g (Gjg, ~ v^m + Q ij n inj ) + g ( $ + * 



t, r— (tq— r)n 



(139) 



Here, for scalar perturbations, 9q*L is given by eq. (134), = — V J u 7 , and — (V'V^ — V 2 ) q, where g, 
negligible during the tight coupling, is a linear combination of 7r 7 and multipole potentials that describe photon 
polarization. The perturbations Ojf^, u 7 , and q are evaluated along the line of sight r(r) = (ro — T)n, assuming the 
observer is positioned at the origin. The "integral visibility function" g(r) is the probability for a CMB photon to 
propagate unscattered from the time r to the present time tq: 



g(r) = cxp 



dr' 



(140) 



Expanding the expression under the line of sight integral (139) over Fourier harmonics in a flat universe, we have 



6(n) 



dr 



d 3 k 

(2tt) 3 



Cin(k) T b (t, k)e 



z(tq — t) n-k 



(141) 



The transfer function Tq in the above equation is constructed from the perturbation Fourier modes, normalized by 
C(tw 0,k) = 1, as 



Tq(t, k) = g 



o 7 + u i — + c l 



d_ 

dr 



d 2 



1 



dr 2 3 



Given that for primordial fluctuations (Cin(k) Cin(k')) = (27t) 3 ^d (k — k') P^(k) and that a plane harmonic projects 



(142) 



onto a spherical one as J d 2 fln Yj* m (h) e m x = 47ryj* Ti (x) i l ji(x), the CMB temperature auto-correlation function C i 
(|6im| 2 ) becomes 



TT 



CT T 



k 2 dk P c (k) 



dr Tq(t, k) ji(k(r - r)) 



(143) 



When CMB polarization or other cosmological 
anisotropics are accessible, additional two-point cor- 
relations can be considered, such as Cf E or Cf E for 
the linear polarization component E generated by 
scalar perturbations, Rcf. [83]. The observed CMB 
polarization can be expressed similarly to eq. (139) 
as a line of sight integral [50] over the perturbations 
that source the photon polarization. Likewise, the 
corresponding Cz's are given by a dk integral of the 
product of the power spectrum P((k) and two time 
convolutions of perturbation variables with ji(k(ro — r)) 
or its derivatives. 



Any contribution to the correlation functions Ci from 
the time of last scattering is characterized by r <C tq. 
The corresponding Bcsscl functions vanish exponentially 
when their argument fc(io — t) ~ fcro is less than I, given 
I 3> 1. Thus the related to the acoustic oscillations con- 
stituents of Ci are essentially affected by only the modes 
that enter the horizon well before the radiation-matter 
equality if I 3> t /S(t e ) « 230, with r e given by eq. (49) 
and WMAP best fit parameters [46]. 
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Neutrino signatures and degeneracies 

For the modes that enter the horizon in the radiation 
era, both the WKB solution (135-136) and the radiation 
era solution (121) should be valid over a positive time in- 
terval 1/k <C r <C r e . Then the comparison of eq. (135), 
where Rf, and x 2 s are negligible in the radiation era, with 
formulas (134, 121) shows that in eq. (135): 

a) The phase shift 6<p is given by eq. (122); 

b) The integration constant A equals 

A = Cin(l + A 7 ). (144) 

The Fourier modes of plasma velocity potential (137) 
acquire the same phase shift and same multiplicative 
change of the amplitude relative to the neutrinoless case. 

The phase shift of the photon acoustic oscillations on 
subhorizon scales can not be produced from adiabatic 
primordial fluctuations by any dynamics involving only 
the photon-baryon plasma and non-relativistic species. 
This is seen easily in real space by noting that the acous- 
tic oscillations correspond to small-scale, appearing sin- 
gular on the Hubble scale, features in the photon Green's 
function. We found in Sec. IV E that in the radiation- 
dominated universe where no perturbations can propa- 
gate faster than the acoustic speed the only such fea- 
tures would be the delta function spikes in eq. (116) at 
x = ±r/V3. 

The spikes would continue propagating away from the 
perturbation origin with the speed c s (t) of eq. (131) past 
the radiation era until the recombination, as, by the gen- 
eral equivalence principle, their local acoustic dynamics 
could not be altered by the gravity of perturbations of 
other species. By the time of recombination, the singular 
part of the photon density Green's function would have 
the form 

d 7 , si „ g (r,x) = D(t) MM - S(t)) . (145) 

The amplitude D(t) depends on the expansion rate of 
the photon background. The calculations of Ref. [56] in 
real space show that D(t) oc e~ k Xs / (1+R.b) 1 / 4 ■ The sin- 
gularity <$D(|a;| — S) becomes 2cos(fc»S) in Fourier space, 
i.e., without neutrinos we would recover the oscillating 
part of the WKB solution (135) with Sip = 0. 

Although the observed period of Ci oscillations de- 
pends on a number of cosmological parameters, such as 
ujb = Qbh 2 , u> m = Qmh 2 , or the total fi, the oscilla- 
tion phase, ideally, can be extracted independently of 
the period. In practice, caution is required. Indeed, even 
in the neutrinoless model the location of the acoustic 
peaks in Ci is not exactly proportional to the wave num- 
ber fc of the extrema of the effective temperature trans- 
fer function (135). This is caused by a variety of effects, 
listed in Sec. 8.3.2 of Ref. [84]. Most of them vanish in 
the I — > oo limit but only as negative powers of I and are 
still sufficiently important to modify quantitative predic- 
tions for Ci phase even at I ~ 3000. Two of the effects 



remain finite for arbitrarily high I. First, the Bcsscl func- 
tion in eq. (143), or a similar equation for polarization, is 
exponentially cut off for I > kr, where r = tq — T 7 dcc 
in the flat model. But (|jz| 2 ) decreases slowly, by a 
power law, for I < kr. As a consequence, an extremum 
of the temperature perturbation transfer function (135) 
at a certain k contributes most to the Ci with I some- 
what less than kr. This shifts Ci peaks toward lower Vs: 
^n'thpoak = (nr/S(T 7dcc ))(n - d> n ), with cj) n > 0. For- 
tunately, the corresponding shift </>„ for Cf T approaches 
a model independent constant value 1/8, Ref. [84]. Sec- 
ond, the phase of Ci oscillations is obscured by the rapid 
decrease of Ci magnitude for I > 10 3 due to the Silk 
damping and smoothing of the anisotropies by the non- 
zero width of the last scattering surface. Nevertheless, 
the numerical analysis of Sec. VI shows the robustness 
of the phase under the change of any standard param- 
eter other than the neutrino density. The phase shift 
signature enables one to constraint N v tightly when a 
sufficiently large I interval is accessible. 

The rescaling of the photon-baryon oscillation ampli- 
tude by 1 + A 7 < 1 for k 1/S(r e ) causes the same 
rescaling of all the photon phase space density and polar- 
ization multipoles, which develop during the decoupling. 
The non-oscillating part of the gravitational potentials, 
also affecting the CMB anisotropy by eq. (135) and (139), 
is however generated by CDM density perturbations and 
rescalcs differently, by 1 + A c > 1, Sec. IV F. In the 
square of the transfer function in eq. (143) some terms 
oscillate in k with the period Afc = ir/S. They come 
from and only from the product of two oscillating pertur- 
bations and must involve the factor (1 + A 7 ) 2 < 1. There 
are terms that oscillate with the period Afc = 2tt/S. 
They are produced by the cross-product of the oscillating 
photon and non-oscillating CDM generated contributions 
and thus get multiplied by (1 + A 7 )(l + A c ), which by 
eqs. (122, 129) is very close to 1. It is hard to find a sim- 
ilar factor for the sum of non-oscillating terms without 
a detailed calculation. These terms consist of both the 
CDM-CDM contributions, scaling as (1 + A c ) 2 , and of 
the non-oscillating parts of 7-7 terms from the squares 
like (cosiys) 2 = i + i cos 2ip, scaling as (1 + A 7 ) 2 . 

Thus we expect that for I ^> r/S(r e ) « 230 the acous- 
tic oscillations of C;'s are decreased by the gravity of 
neutrino perturbations by the factor (1 + A 7 ) 2 . But the 
magnitude of the Ci oscillations depends on other phys- 
ical quantities as well. 

First of all, on the primordial power spectrum P^(k). 
In principle, this unknown could be excluded by compar- 
ing the CMB anisotropies with the matter density per- 
turbations. The latter scale differently under N v vari- 
ation, as seen in Sec. IV F and in the next subsection. 
However, given the large experimental error on the mat- 
ter perturbations, a better method to exclude Pq is to 
compare the height of the initial acoustic peaks, enter- 
ing the horizon closer to the matter era and hence less 
affected by neutrino perturbations, to the height of the 
subsequent peaks. This signature is degenerate with the 
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power spectral index n s = dhiP^/(dhik) + 4 and its run- 
ning in k, especially for only a limited number of observed 
peaks [34]. 

Second, the scaling (1 + A 7 ) 2 is directly applicable 
only if one compares the models in which the photon 
subhorizon dynamics is identical in angular and rcdshift 
coordinates 12 since the radiation era until the present. 
Here is the list of the dimensionless quantities that char- 
acterize the cosmological background expansion and the 
local photon-baryon dynamics: 13 1. The ratio of the 
universe expansion rates at any redshift z and at the 
present: H(z)/Hq. It shows how distances are mapped 
to angles: (z,k/h) — > I. It also affects the photon and 
baryon perturbation dynamics, considered in the dimen- 
sionless variables dz = —H(z) dr and Hodx. This ratio, 

H(z)/H = [il p(z)/p(0) + (1 - fi)/a 2 ] 1/2 , (146) 

is specified by today's total density parameter ~ 
f2 m + f2dc, Cl m , the redshift of the radiation-matter equal- 
ity 1 + z oq = p m fi/ Prfiy and the dark energy "equation 
of state" Wde(z) = pdc/pdo- 2. The ratio of baryon 
and photon densities. Pb/p-y, controlling the photon- 
baryon plasma dynamics. 3. The ratio of the photon 
free- flight time r c (r), eq. (132), to the particle horizon 
size r. This ratio determines the integral visibility func- 
tion g(z), eq. (140), and the important Silk damping 
scale, eq. (136). 

The present photon density p 7i o is well constrained by 
the COBE measurement of CMB temperature 2.725 ± 
0.002 K, Ref. [74]. It is expected to redshift predictably 
as pj = (1 + z) p-yfi deep into the radiation era. The to- 
tal radiation density depends on N v as p r = p 7 (l + a v ), 
where a v ~ 0.23 iV„, eq. (44). Given this, two mod- 
els with different N v will have same ratios pb/p-y and 
H(z)/Ho, and the related p m ,o/Prfi = 2 oq + 1, if these 
models have the same ujb = ilbh 2 , CI, fi m , and Wd e (z), 
but their Hubble constants scale as h cx \/l + a v . 

The two top panels in Fig. 5 show the relative change 
SCi/Ci for CMB temperature and polarization spectra 
under N v variation from 2.5 to 3.5. The solid curves cor- 
respond to preserving the parameters u>b, fl = 1, f2 m , 
Wdc = — 1, and h/y/1 + a v , as described above, as well 
as the primordial power spectrum and the primordial he- 
lium fraction Y = pn e / Pb- As seen from the plots, the 
model with a larger N v has noticeably stronger damping 
on small scales. Indeed, for a fixed ionization fraction 
x e = n e /riB Tb e Ju>H, where ujr = (1 — Y)ujb specifics 



A universe with larger neutrino density expands faster, at least, 
in the radiation era. The corresponding photon temperature 
must, therefore, decrease faster in the time coordinate. 
We yet ignore the somewhat relevant to CMB ratio pj,/p m . As 
seen in the next subsection, this ratio is very important for mat- 
ter evolution. But it affects CMB anisotropics rather mildly, 
through the non-oscillating CDM potential term in eq. (135) or 
through the CMB lensing by matter structure. 



the primordial hydrogen density, the ratio 
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increases with N v as h cx y/\ + a v . 

The full conformal equivalence of the photon subhori- 
zon dynamics can be straightforwardly implemented in 
numerical calculations. The dotted curves on the top 
panels of Fig. 5 show SCi/Ci when the global cosmolog- 
ical parameters vary as before but the numerically inte- 
grated equations in both models use the same r c /r and 
the same integral visibility function g(z). The remaining 
Ci change for the modes I 3> 200, entering the horizon 
before the radiation-matter equality, is uniform power 
suppression and a constant phase shift, as predicted for 
the effects of neutrino perturbations. The shift I — > I + SI 
in the numerical calculations is only 80% of the expected 
51 = Alp Ca bSLp/ir ~ 4.6, for the period of the acoustic 
oscillations A/ pca k ~ 300 from Ref. [85] and Sip evalu- 
ated from the radiation era result eq. (122). Likewise, 
the power suppression is somewhat less than the radia- 
tion era prediction 5(1 + A 7 ) 2 ~ —4.8%, with A 7 given 
by the leading term in eq. (122). These discrepancies are 
caused by the residual effect of the large-scale relativistic 
correction, c.f. the first term in eq. (117), non-negligible 
matter density during the horizon entry, and 0(i? 2 ) cor- 
rections in eq. (122). 

As suggested by eq. (147), the ratio t c /t in the com- 
pared models could be matched by varying the helium 
abundance Y. Apparently, when all the helium has re- 
combined but the hydrogen remains fully ionized (x e = 
1), the quantity (147) remains constant if Y varies as 
(1 — Y) cx h cx VI + a v . It is less obvious that the after 
such a rescaling, t c /t coincides in both models during the 
recombination, when x e violently changes with time. De- 
spite the complicated nature of recombination, the phys- 
ical mechanism that is primarily responsible for x e evo- 
lution around the peak of the photon visibility function g 
does lead to degenerate x e (z) and so degenerate t c /t. 

At the redshifts 1400 > z > 800, which are of the 
most interest, the dominating process leading to H re- 
combination to the ground state Is is the suppressed 
2s — > Is + 7 + 7 decay. Faster transitions to Is with the 
emission of a single photon do not create more hydrogen 
in the ground state because the emitted resonance pho- 
ton soon ionizes or excites another H atom in the ground 
state. Using the approximations of Ref. [86] (see Ref. [87] 
for a recent review) , one can find the rate of x e change by 
assuming that due to the overcooling of the universe from 
the delay in recombination the fraction of all the excited 
H atoms is negligible comparatively to the concentration 
of H in the ground state x\ s = n\ s jnn and to x e . Then 
x\ s + x e ~ 1, hence dx\ s ~ —dx e . The change in x\ s is 
mainly due to the spontaneous decay 2s — > ls + 7 + 7 
with the lifetime A ~ 8.2s _1 : dxi s /dt ~ Ax2 S - (The 
photo-excitation Is — ► 2s is negligible when z < 1300.) 
x e decreases in ep recombination and increases in hy- 
drogen photo-ionization. With realistic approximations, 
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FIG. 5: The relative change in Cj T (top left), Cf B (top right), matter density perturbation 8 m (k/h) (bottom left), and Cf K 
(bottom right) when N v varies from 2.5 to 3.5. The solid curves show the changes when all the other parameters, listed in 
Sec. VI, are fixed. The top two panels also show the change for fixed recombination history and equivalent Silk damping 
(dotted, blue). The dotted, blue curve on the bottom left panel gives the change in 5 m when ujb/co m is held fixed. 



Ref. [86], one finds: dx e /dt = — ax e x p nn + (3x2 S , where 
x p ~ x e , and a and /3 at a given CMB temperature, de- 
termined by the redshift, are independent of cosmological 
parameters. The elimination of X2 S from the equality of 
dxi s /dt and —dx e /dt yields: 

dx e aA 

It ^ ~ JTI XeXpnH ■ (148) 

When both d/dt = —(1+z) H(z) d/dz oc h and nn oc (1 — 
Y) uji, vary in proportion to ^/l 4- a„, as discussed above, 
eq. (148) predicts that the function x e {z) is unchanged, 
and so is unchanged the last ratio in eq. (147). 

Numerical calculations show that this conclusion holds 
very well after the end of helium recombination at z ~ 



1500. For a variation of N v from 2.5 to 3.5 and the cor- 
responding adjustment of Y by —0.051, the ionization 
fraction x e changes at most by 0.5%, and even twice as 
little within the peak of the photon visibility function g. 
(For comparison, without the adjustment of helium abun- 
dance, the change of x e reaches 6% within the visibility 
peak.) At 2000 < z < 5000, when helium is singly ion- 
ized, the decrease of the helium density decreases the 
total number of free electrons by about 2%, and more 



Astrophysical constraints on the primordial helium abundance 
vary among groups and span a range Y = 0.238 ± 0.010 (2cr), 
Refs. [6—10], a review in Ref. [88]. Taking this prior, Y needs to 
be considered as a free parameter in a CMB analysis only for the 
accuracies lAAVI < 0.4. 
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at higher redshifts, when the second helium electron un- 
bounds. But the related increase of the Silk damping 
scale, eqs. (132,136), is insignificant, and the effect is 
quite negligible for dg/dz. The change SCi/Ci for TT 
and EE spectra under the considered variation is seen 
from Fig. 6 (solid line), ft is compared with the the- 
oretical reconstruction (dashed line) obtained by shift- 
ing Ci by 51 = 0.8AZ pca k<5<y9/7r ~ 3.7 and rescaling it by 
5(1 + A 7 ) 2 = -4.8%. The plots show that the Y adjust- 
ment removes most of the damping seen in Fig. 5. 

The approximate conformal degeneracy of CMB dy- 
namics among the models with changed dynamical time 
and length scales should be distinguished from the well 
known degeneracy of primary CMB anisotropies under 
any variation of H(z) for z < z 7 dcc that preserves the 
angular diameter distance to the last scattering surface, 
Ref. [59]. The latter degeneracy is generally violated by 
the late ISW effect. The found one is respected by ISW 
but is somewhat violated by the gravity of matter per- 
turbations, c.f. footnote 13 on p. 22, and, of course, by 
the gravity of neutrino perturbations. 

The internal CMB dynamics is shown to be confor- 
mally invariant in an increment of N v if the helium abun- 
dance Y is decreased as AY/AN V ~ —0.05. In contrast, 
the standard Big Bang nucleosynthesis (BBN) then pre- 
dicts increased 4 He production. 15 Indeed, for AN V > 
the universe expands faster and less neutrons decay by 
the time it cools sufficiently to allow their conversion into 
helium. Of course, one should not presume a priori that 
the ratio of neutrino and photon energy densities is pre- 
served since the nucleosynthesis until the matter era, i.e. 
that N v in BBN and CMB physics are equal. But if they 
are not, the non-degeneracy of CMB dynamics along the 
BBN-predicted curve Y(N V ) makes it easier to spot the 
discrepancy. 

Neutrino perturbations can affect the gravitational po- 
tentials even after the photon decoupling if radiation is 
still a significant component of the total energy density. 
When the gravitational potentials are time dependent, 
there is a contribution to the CMB anisotropy (139) that 
depends on the line of sight integral of <f> + ^, the so- 
called integrated Sachs- Wolfe (ISW) effect. However, as 
discussed in Sec. Ill C, for superhorizon modes the deriva- 
tive 9($ + ^)/dR v changes insignificantly after recom- 
bination in the standard models with matter radiation 
equality at z ~ 3000. Thus the modification of the ISW 
effect by neutrino anisotropic stress does not play a ma- 
jor role as a neutrino signature in CMB. Specifically, on 



15 The sensitivity of the BBN yield of 4 He and D/H to variation 
of N v and jj = n B /n 7 about N v = 3 and the WMAP [13] "best 



fit" value r\ ?s 6.1 X 10 is roughly 
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large scales (I < ?xst pcak/3 — 70) the ISW effect becomes 
insensitive to perturbations of relativistic neutrinos. 



B. Matter power spectrum 

The growth of the acoustic perturbations in the tightly 
coupled photon-baryon fluid is halted by photon pressure, 
described in eq. (7) by the term — c 2 V 2 <i 7 , where c s > 1/2 
before recombination. Likewise, relativistic neutrinos 
and CMB photons after their decoupling are stabilized 
against gravitational collapse by the effective pressure, 
described by the same term with c a = 1/ -\/3, due to the 
velocity dispersion. In addition, as noted in the intro- 
duction, the perturbation modes of free-streaming parti- 
cles on subhorizon scales decay by "directional" damping. 
Assuming the dark energy too does not strongly clus- 
ter on small scales, only CDM and after recombination 
baryons can cluster sufficiently to generate non-negligible 
gravitational potential inside the subhorizon. On sub- 
horizon scales eq. (23) simplifies to V 2, 3> = •yd. Then 
eqs. (7, 24) applied to non-relativistic species n with neg- 
ligible pressure (c 2 ~ 0) give 



d n + Hd n = 



777 



(subhor) . (149) 



One can distinguish the following stages of the linear 
matter perturbation growth from their entry of the hori- 
zon in the radiation era to the present: 

1. Growth of CDM perturbations in the radiation- 
matter universe until the recombination while 
baryons are tightly coupled to CMB. 

2. Decoupling of the baryons from CMB and joining 
the matter gravitational collapse at Zd ~ 1090. 

3. Growth of the pressureless CDM-baryon matter 
perturbations through the subsequent universe evo- 
lution, affected by dark energy and, possibly, global 
curvature. 

Eq. (128) for CDM density perturbations in the radi- 
ation era on subhorizon scales can be presented as 



d ( J' X<<T) (T,k) = 
= -6Cin(l 



(150) 



A c ) <A + In 



(1 + Sc)- 



where A is the same function of rh for all k and R u : 
rh\ 1 



A(rh) = In 



(rad. era) . (151) 



At later times the matter density perturbation must be 
of the form 



dm{T,k) =-6Cin(l + A c )<M + Bln 



(1 + 5c) 



,(152) 



(based on the numerical results from Refs. [11, 89, 90].) 



where A and B depend on rh and on the cosmolog- 
ical parameters affecting the background evolution at 
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FIG. 6: The relative change in Cf T (left) and C; BB (right) for N v variation from 2.5 to 3.5 and the adjustment of the helium 
abundance Y by —0.051 (solid). It is compared with the change from rescaling Ci's by —4.8% and shifting them to smaller l's 
by 51 = 3.7 (dashed). 



low redshift. We argue that throughout linear evolution 
A and B are independent of the mode wave-number k 
and of i?„, provided that the compared models agree 
in f2 m , in Of, , and in the background effective equa- 
tion of state w(z) = p/p. This would immediately fol- 
low from the linearity of the evolution with the initial 
conditions (150) in the radiation era if the equations 
of matter dynamics in terms of the dimcnsionless vari- 
able tHq oc rh coincide. The latter is manifestly so by 
cqs. (149,42) for the stages 1 and 3 of the three stages 
of linear matter clustering identified above. As for the 
stage 2, the scale and R v invariances of the dynamics are 
violated by the small oscillations of the baryon density 
and velocity at the decoupling due to the acoustic oscil- 
lations in CMB, Ref. [2]. However, this effect becomes 
negligible at sufficiently small scales because of the Silk 
damping. Hence we assume that the baryons decouple 
with initially negligible density and velocity perturba- 
tions. Then all the stages of the matter evolution are 
scale- and ^-invariant. Then, from eq. (152), the mat- 
ter density perturbation modes d m {z,k, R v ) can be ob- 
tained from the modes in the model with R v — > 16 and 
the identical 17^, and w(z) as 

d m (^z,^,R u ^=(l + A c )d£»^z,(l + 5c) ] ^,(153) 

where for the considered small scales, entering the hori- 
zon in the radiation era, A c and Sc are given by eq. (129). 

This effect of neutrino perturbations on the matter 
density (dotted curve on the left bottom panel of Fig. 5) 
is too small to provide by its own a useful informa- 
tion about N v from the available cosmological probes. 



16 In the models with dark matter, cosmological constant, global 
curvature, but negligible baryon density (pb/pc -C 1), the linear 
matter density perturbation well after the radiation era on the 
small scales is 
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where d nl = — 3£i n . Analytical, albeit more complicated, expres- 
sions also exist for non-negligible baryon density [2, 82] . 



A realistic data analysis constraining the abundance 
of ultra-relativistic neutrinos should also include CMB 
data. However, the variations of CMB spectra with N v , 
Sec. VA, are less contaminated by neutrino unrelated 
physics if u>b = 0/,/i 2 , rather than fib, is fixed. Likewise, 
the Big Bang Nucleosynthesis (BBN) constraints uib and 
not Of, . Thus it is more practical to consider the vari- 
ation of matter density modes with in the direction 
of the maximal CMB and BBN degeneracies - under a 
fixed ujb. 

The magnitude of matter perturbations after the radi- 
ation era is sensitive to the parameter 
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e.g. Ref. [2], which is currently estimated as 0.17 ± 
0.01 [13]. One reason for this sensitivity is the slower 
growth of CDM perturbations prior to the recombina- 
tion in the model in which a larger fraction (3 of the non- 
relativistic matter is withheld from gravitational collapse 
by the photon pressure. The other reason is a greater di- 
lution of the growing CDM perturbations by almost un- 
perturbed baryons when the latter finally decouple from 
photons. 

It is apparently impossible to change the neutrino den- 
sity preserving all of the parameters u> m /ui r , u>t>, and (3. 
If one varies R v while keeping 1 + z cq = w m /oj r and Wb 
fixed, to minimize the changes in CMB spectra, then the 
parameter (3 will vary. For example, (3 decreases by ap- 
proximately 41% for N v change from to 3, and by 14% 
for N v variation by 1 around its standard value 3.04. 
This is yet another source of breaking the degeneracies 
between N v and the other cosmological parameters. The 
significance of (3 variation for the growth of matter per- 
turbations is evident from comparing the dotted (fixed Qb 
and Q m ) and solid (fixed Ub and f2 m ) curves in Fig. 5. 



VI. FORECASTS FOR FUTURE 
EXPERIMENTS 

In this section we use numerical solutions from CMB- 
FAST and apply them to predict the precision to which 
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N v can be constrained from future experiments. We 
follow the basic approach of Ref. [91]: we evaluate the 
standard error in a cosmological parameter s; as Asj = 
(cc -1 )*/ 2 , where a. is the Fisher matrix: 

«-EE? c --w,cr)f • (155) 

I X,Y 1 J 

Here, Cov -1 is the inverse of the covariance matrix, Sj 
are the unconstrained cosmological parameters, and X, Y 
stand for the observable power spectra. We limit the 
analysis to CMB spectra TT, EE, TE, and lensing con- 
vergence power spectrum kk as measured from CMB it- 
self. For each I, one has to invert the covariance matrix 
and sum over X and Y. The derivatives were calculated 
by finite differences. The step was usually taken to be 
about 5% of the value of each parameter and symmetric 
around the pivot point at the best fitted WMAP model, 
Ref. [13]. 

An experiment is fully characterized by its sky cover- 
age /sky, temperature, polarization, and convergence in- 
strument (or reconstruction) noise w^, , Wp, and w" 1 
and by the beam smearing window function B7 = 
e z(i+i)0£/(8in2)^ yrffa Q b measuring the width of the beam. 
For example, the covariance element for TT is given by 

Cov (dr, CD = * , , {CJ T + W) a - (156) 

The full covariance matrix for the CMB power spectra 
is given in Ref. [92]. For the lensing convergence noise 
spectrum we take the values from the maximum like- 
lihood method developed in Ref. [93], which gives the 
lowest reconstruction noise. In addition, we impose a 
maximum Z max = 3000 cutoff on all the spectra and we 
do not include the information from higher I. The jus- 
tification for this is that scattering off moving electrons 
during and after reionization leads to an additional com- 
ponent in the CMB that cannot be separated from the 
CMB on the basis of frequency information. This com- 
ponent is rather uncertain since it receives contributions 
from perturbative and nonlinear structures, as well as 
from the patchiness in ionizing fraction. It is expected 
to be significantly less important for polarization, so our 
approach may be conservative for polarization sensitive 
experiments. 

In our analysis we use parameters for several future ex- 
periments. While WMAP constraints on N v arc not yet 
competitive with the nucleosynthesis limits, Planck satel- 
lite will improve the sensitivity of WMAP by an order of 
magnitude. In parallel, there will be ground based exper- 
iments, such as SPT 17 or the proposed ACT 18 that will 
extend Planck to smaller angular scales. Whether or not 
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they will be polarization sensitive remains to be deter- 
mined, so we explore both possibilities. Finally, we also 
explore the prospects of an ambitious high resolution, 
low noise satellite dedicated to polarization (CMBPOL), 
which will be able to measure with high accuracy not 
only CMB temperature and polarization but also matter 
power spectrum using the weak lensing induced signal in 
higher order correlations. We also explore if adding in- 
formation from more local probes, such as weak lensing, 
can further reduce the uncertainties. 

We explore 11 parameters in our analysis. These are 

massive ; 

P,n s ,n' s ,Y}. 

The first one is proportional to the ratio of matter 
density to radiation density (photons plus neutrinos), 
the second to proper baryon density, and the third to 
the proper density of massive neutrinos. Side is the dark 
energy density relative to the critical density, w dc is 
its effective constant equation of state Pde/pde, T re ion is 
the reionization optical depth, P is the amplitude of 
curvature power spectrum at k = 0.05/Mpc, and n s and 
n' s = dn s /d\nk are the scalar spectral index and its 
running at k = 0.05/Mpc. We do the analysis for both 
fixed and unconstrained helium abundance Y. 

We ignore tensor perturbations since their contribution 
is limited to large scales, where neutrinos do not play a 
major role. We also ignore BB power spectrum, which is 
useful primarily as a tracer of tensors. (On small scales it 
can be a useful tracer of matter power spectrum and can 
provide additional constraints on the convergence power 
spectrum as extracted from the nongaussianitics in CMB; 
we ignore this additional information here.) We use the 
basis functions dCf jdsi as defined by these parameters, 
but we also explore any possible numerical instabilities 
due to the variation in the angular scale of the acoustic 
horizon, which is affected by variations in several of the 
parameters. 

Fig. 5 shows the N v derivatives of Cf for X = TT, 
EE, kk, and of the matter fluctuation S m . The deriva- 
tives of CMB spectra with lensing arc plotted, but the un- 
lensed ones are qualitatively similar and the wiggles are 
not due to lensing, but due to the phase shift (the lensed 
derivatives were used in the Fisher matrix analysis). The 
solid curves on the C; plots show the N v derivatives with 
the other parameters kept fixed. As discussed in previ- 
ous sections, changing the effective number of neutrinos 
changes not only the initial phase and amplitude of the 
acoustic oscillations but also the angular scales of the 
acoustic horizon and oscillation damping. To separate 
these effects, we also show with the dotted, blue line the 
change in the spectra while keeping the visibility func- 
tion and the photon free flight length in the units of H^ 1 
unchanged at every redshift. 

The modification of CMB spectra due to the change in 
the angular scale of acoustic oscillations is described by 
multiplicative rescaling I — > (1 — 0.002 AN u )l. By itself, it 
is degenerate with other effects that change the angular 
acoustic scale, such as variation of Wdc- The combined ef- 
fect of the additive and multiplicative phase shifts is such 



27 




1000 



aooo 



3000 




3000 



FIG. 7: The derivatives of Cf T (top) and Cf E (bottom) with 
respect to w& (solid, red), u m /uj r (dotted, black), dn a /d\nk 
(short-dashed, cyan) and n a (long dashed, blue) and Y 
(dotted-dashed, green). All were shifted to match the angu- 
lar scale of acoustic horizon. This explains lack of oscillations 
when varying uo m /u r or Y . Varying uJb changes the ampli- 
tude of the oscillations, so adjusting the angular scale does 
not make their change vanish. Note that the variation in Y 
appears as a change in damping length. 



that the phase shifts cancel exactly at I ~ 1500, but not 
at other values of I. For the temperature anisotropics, 
the phase shift is barely visible since temperature os- 
cillations are weaker due to the competing effects from 
density and velocity terms. They are further suppressed 
by lensing. For polarization, which has more prominent 
acoustic oscillations, the phase shift remains visible and 
can be clearly distinguished from the change in angular 
size of the acoustic horizon. This suggests that polariza- 
tion information is crucial in extracting neutrino signa- 
ture. Quantitative analysis in Table III confirms that. 
Note that 3Ci/dN v approach zero at low I. There is no 
significant neutrino dependent contribution coming from 
the integrated Sachs- Wolfe term at the low /, as discussed 
in Sees. IIIC and VA. 

In Fig. 7 we show the derivatives of Cj T and Cf E with 
respect to some of the other parameters. The displayed 
derivative are for u> m /u> r , uib, n s , dn s / 'dhik, and Y and 
they are taken by keeping the angular scale of the acous- 
tic horizon constant. Other derivatives are trivial and 
we do not show them here: the amplitude and reioniza- 
tion optical depth only rescale the amplitude, while dark 
energy density and its equation of state only change the 
angular scale of the acoustic horizon and make no ef- 
fect after this change is taken out. Massive neutrinos 
have a minor effect on CMB over the range of interest. 



It is evident that none of these can mimic the additive 
phase shift generated by neutrinos. This signature thus 
uniquely identifies the presence of neutrinos in the con- 
text of adiabatic initial conditions. 

If matter power spectrum information is available, one 
can exploit the fact that while neutrino perturbations 
suppress CMB anisotropies, they enhance the matter 
power spectrum, Figs. 4 and 5. For AN V = 1 at a fixed 
UJb, k/h = IMpc -1 , and z = this gives AP m /P m w 
0.12, of which 0.02 is due to neutrinos alone and the rest 
to the variation of tOb/uj m , as discussed in Sec. VB. At 
the same time, the CMB power spectra are suppressed 
by 15% at I = 3000 if Y is fixed. Whether or not this is a 
useful method to break the degeneracies depends on the 
accuracy with which matter power spectrum can be ex- 
tracted. We performed the analysis for the experiments 
in Table III assuming weak lensing information as can be 
extracted from the CMB itself. A Fisher matrix analysis 
with lensing reconstructed convergence power spectrum 
using minimum variance maximum likelihood errors [93] 
fails to improve significantly the accuracy on N v , even 
with the polarization information, which allows for a bet- 
ter reconstruction of the convergence spectrum. This is 
in contrast to other parameters such as massive neutri- 
nos, for which lensing information significantly improves 
the bounds [94]. 

Unless the astrophysical constraints on the primordial 
helium abundance Y will improve by an order of magni- 
tude, Y should be included in the list of unknown cos- 
mological parameters for the accuracies of the order of 
AN V ~ 10 _1 . On the other hand, if we assume that the 
neutrino fraction of the radiation is unchanged between 
BBN and recombination and there are no significant v jv 
asymmetries then BBN tightly limits the helium frac- 
tion, to AF ~ 10~ 3 or better [11, 12]. For this reason, 
we performed the parameter prognosis with and without 
Y parameter. As found in Sec. V A, an adjustment of Y 
can compensate the changes in the acoustic and damping 
scales due to neutrino density but it preserves the phase 
shift signature due to the gravity of neutrino perturba- 
tions. Correspondingly, the CMB limits on N v somewhat 
broaden but remain tight even with the unknown value 
of Y. Interestingly, the CMB data itself can be used 
to constraint the primordial helium abundance indepen- 
dently of the astrophysical measurements and the related 
systematic uncertainties. For the considered CMBPOL 
experiment with the polarization data included we find 
a(Y) ~ 0.005 from the CMB alone. 



VII. CONCLUSIONS 

In this paper we study analytically the evolution 
of cosmological perturbations in the presence of ultra- 
relativistic neutrinos. While dynamical equations for cos- 
mological perturbations have been known for a while [3, 
57, 95-97], their analytical solutions exist only in a hand- 
ful of cases and are restricted to the fluid description. 
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TABLE III: Standard deviations on N„ as expected from Planck, ACT, Planck+ACT and CMBPOL using temperature data 
only (TT) and added polarization (TT + TE + EE) . The primordial helium abundance Y is considered a priori unconstrained 
for the last column results and fixed by independent measurements in the preceding two columns. Adding weak lensing 
convergence as reconstructed from CMB (TT + TE + EE + kk) does not significantly improve the bounds, even assuming 
polarization information is available. 



The best known examples, e.g. Rcf. [60], are the solu- 
tions for CDM and photon-baryon plasma in the matter 
and radiation eras or in the subhorizon limit, and for 
supcrhorizon metric perturbations. In contrast, neutri- 
nos cannot be modeled by a fluid and their phase space 
distribution should be considered. 

Most of the recent publications abandoned the ana- 
lytic approaches and relied on numerical results from 
Boltzmann integrator codes. While, in principle, there is 
nothing wrong with this, analytic solutions often lead to 
deeper understanding of the problem that can reveal the 
new directions of exploration. They sharpen the focus 
on the features that are unique and cannot be mimicked 
by the variation of other parameters. Care must be exer- 
cised when performing numerical analysis and parameter 
forecasting for future experiments. The computational 
errors must be well controlled, otherwise they can lead 
to artificial breaking of degeneracies. Besides, the param- 
eter space of forecasting is often small and with the ad- 
dition of new parameters new degeneracies may open up. 
For example, while so far only simple parameterization 
of the primordial power spectrum have been explored, 
one could consider its more general parameterization, in- 
cluding the running of the slope, its running, etc. as free 
parameters. In this case analytic solutions can provide 
a better understanding of whether the limits on a given 
parameter are robust against adding new parameters. 

We obtain analytic perturbativc in p v jp r solutions 
for cosmological perturbations in the presence of ultra- 
relativistic neutrinos. Much of the success is due the 
following result, equally useful for fluid models. We find 
that a simple redefinition of the independent dynami- 
cal variables that is consistent with their classical in- 
terpretation and preserves them on small scales, elimi- 
nates all the time derivatives of the non-dynamical metric 
perturbations from the evolution equations in the New- 
tonian gauge. The resulting description of cosmologi- 
cal perturbations acquires the advantages of an initial 
value Cauchy's problem while remains formulated in the 
Newtonian gauge, which is fully fixed and is especially 
suitable for describing the physics of CMB. Moreover, it 
turns out that even in the solvable fluid models the solu- 
tions for matter or radiation density perturbations, c.f. 



eqs. (117, 127), appear far simpler in the redefined vari- 
ables. In addition, these variables are generally constant 
on supcrhorizon scales. 

While most of the previous literature has focused on 
Fourier space analysis, we also consider perturbation 
Green's functions in real space [55, 56]. The latter be- 
come indispensable for the analytical study of neutrino 
perturbations. They also allow one to prove quickly that 
without neutrinos the cosine form of the acoustic density 
oscillations in the radiation era cannot be modified by 
the gravitational feedback processes. 

We use the zero, in the powers of p v j p r , order solutions 
for neutrino perturbations to derive the analytic expres- 
sions for the CMB and matter density fluctuations in the 
linear order. We show that these first order solutions 
are for the most part sufficient for a quantitative inter- 
pretation of numerical solutions. Finally, we use the full 
numerical solutions from CMBFAST to derive parameter 
forecasts for various planned experiments. The presented 
methods can be straightforwardly extended to other ap- 
plications such as tensor modes and massive neutrinos. 
We plan to address some of these in future work. 

The distinctive cosmological features of ultra- 
relativistic neutrinos are due to their free streaming at 
the speed of light. The free streaming creates neutrino 
anisotropic stress, perturbing Newtonian metric even for 
supcrhorizon modes. In real space, it also leads to per- 
turbing the photon-baryon plasma beyond the acoustic 
horizon of a primordial perturbation. In Fourier space 
this manifests itself as the phase shift in the acoustic 
oscillations that is generated at horizon crossing. This 
phase shift is unique in the sense that for adiabatic per- 
turbations no non-relativistic or fluid matter can gener- 
ate it. The effect changes the phase additively and cor- 
responds to ~ —4 for AN V = 1. In contrast, any 
change in the angular scale of acoustic horizon acts as 
multiplicative rescaling I — » cd. The two shifts are only 
degenerate at a single I and can be distinguished in gen- 
eral, Fig. 5. The effect is more visible in polarization, 
which has sharper acoustic peaks relative to temperature 
anisotropy, where the density and velocity contributions 
to the Ci oscillations partially cancel. As a result, the 
precision of determining the effective number of neutrino 
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species can be improved dramatically if polarization in- 
formation is included. 

Phase shift is not the only signature of neutrinos in 
CMB. The free streaming neutrinos also suppress the os- 
cillation amplitude of the CMB modes entering the hori- 
zon in the radiation era. A change ANl, = 1 leads to 
ACi/Ci ps -0.04. Since the CMB modes entering in 
the matter era and not experiencing the suppression are 
limited to large scales, where sampling variance is large, 
this effect by itself cannot be extracted with high preci- 
sion. However, neutrino perturbations amplify the CDM 
modes entering the horizon in the radiation era. The 
effect is further enhanced by the fact that while CMB 
physics is more sensitive to the ratios LUb/uj 7 and uj m /uj r , 
specifying the acoustic dynamics and the background 
evolution, matter fluctuations are also sensitive to the 
ratio LOt/tUnn which, if fixing the natural CMB variables, 
cannot be held fixed under varying N u . This changes the 
present matter fluctuation spectrum on small scales by 
AP m /P m PS 0.12 for AN U = 1, of which 0.04 is due to 
neutrinos and the rest to u>b/<x> 7 variation. 

It is unclear how accurately can this effect be extracted 
from local probes of large scale structure, such as galaxy 
clustering and weak lensing, since nonlinear evolution 
will complicate or, in the case of galaxy clustering, pre- 
vent its determination on small scales. Weak lensing 
of CMB traces matter fluctuations on larger scales and 
higher rcdshifts than any other method. It may be the 
optimal tool to use here since nonlinear reconstruction 
methods using the nongaussian information, especially 
from polarization data [93, 98], can achieve high signal 
to noise on the projected convergence power spectrum. 
However, from the Fisher matrix analysis we find that 
lensing of CMB cannot improve the limits from primary 
CMB and its polarization significantly. 

Finally, N v variation changes the relativistic energy 
density and thus changes the relation between the ex- 
pansion factor and time. This leads to a change in the 
proper size of the acoustic horizon and so in its angular 
size, which determines the positions of acoustic peaks. 
The angular size of the horizon, however, is degenerate 
with other parameters, such as those changing the angu- 
lar diameter distance to recombination. The change of 
the expansion time scales also modifies the recombination 
process, the visibility function, and the angular damping 
scale. The effect on the CMB power spectrum can be 
significant, reaching 15% power suppression at I = 3000 
for AA„ = 1, Fig. 5. However, this can be mimicked 
by different primordial helium abundance: a change of 
AN V = 0.1 is compensated by AY ~ -5 x 10" 3 . If 
CMB data is used to constrain ns/n 7 at BBN then the 
standard BBN limits on Y are already at the level of 
AY < 10~ 3 [11, 12], suggesting that Y can be assumed 
fixed. These limits are not applicable in the models where 
the photon entropy changes between BBN and CMB de- 
coupling or in non-standard BBN models with v/v asym- 
metries or particle decays. 



In summary, the effects of ultra-rclativistic neutrinos 
on CMB and matter power spectrum are generally small. 
This is why only weak limits on the neutrino background 
density have been placed from the available observations. 
On the other hand, neutrinos give rise to unique effects 
which exist on small scales and are thus less limited by 
sampling variance. As a consequence, future CMB ex- 
periments should be able to improve the limits signifi- 
cantly. While Planck will be able to determine N v with 
a standard deviation 0.24, or 0.20 if Y is constrained, a 
dedicated CMB polarization experiment should improve 
this bound even further, reaching accuracy levels of 0.09 
without Y constraint, or 0.05 if Y is constrained. This 
will allow one to test the details of neutrino decoupling 
and the scenarios giving rise to a nonstandard number of 
neutrino species. 

Acknowledgments 

SB is grateful to E. Bertschinger, J. R. Bond, and 
A. Loeb for their encouragement in studying neutrino 
perturbations, and to C. Lunardini for a discussion of 
BBN physics. We both thank E. Bertschinger and 
V. Mukhanov for useful discussions, A. Makarov and 
P. McDonald for valuable help with numerical calcula- 
tions, and S. Weinberg for suggestions allowing to extend 
the supcrhorizon conservation laws of Sec. II B to the 
inflationary epoch. SB has been supported by Prince- 
ton University Dicke Fellowship; US by Packard Foun- 
dation, Sloan Foundation, NASA NAG5-1993 and NSF 
CAREER-0132953. 

APPENDIX A: COSMOLOGICAL DYNAMICS 

Given mutually non-interacting, except than gravita- 
tionally, groups of cosmological species {a}, one can de- 
fine 19 for each group an energy-momentum tensor Tg v 
that satisfies the local conservation law T£ u . v = 0. The 
total energy-momentum tensor 

a 

sources space-time curvature perturbations, as described 
by the Einstein equations R^ v — \g^ v R = 8nGT^ u . 

1. Background 

We study cosmological perturbations relative to a spa- 
tially homogeneous and isotropic model with the metric 

ds 2 = a 2 (r) (-dr 2 + 7ydxW) . (A2) 



The definition is Tq V (x) = (2/y/—g) 8S a /^g^v (x), where S a are 
the terms of the action that describe the species a. 
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The spatial part of the background metric may be written 
as 



Ho 



dx t dx 3 



dr 2 



1-Kr 2 



r 2 dn 2 



(A3) 



For most applications of this paper, except for the late 
time evolution of the matter perturbations, the back- 
ground curvature K can be neglected. In this case, we 
take 7ij = 5 i3 . 

The background expansion rate with respect to the 
conformal time r is denoted by Tl = a/ a = aH. By the 
Fricdmann equation, it equals 



H 2 = 



8irGc 



p-K, 



(A4) 



where p = —T°q. 

In the unperturbed universe, 

T a % = - Pa , t^ = o, v;: ; 5) Pa 



(A5) 



Assuming that the species pressure p a is uniquely speci- 
fied locally by the species energy density p a , we introduce 



Pa 

Pa ' 
dp a 
dp a 



dw a 

d In p a 



(A6) 
(A7) 



For the similar quantities applied to all the cosmological 
species together, with p = J2 a P a anc ^ P = J2 a P a ' we 
have 



P 
P 

P 
P 



_ \ ~« x a 

l + W ~ ^ 1 + W a 

a 



Here 



<)}) 

^Pj adiab 



Pa +Pa 
P + P 



(A8) 
(A9) 

(A10) 



are species enthalpy abundances, satisfying J2 a x a = 1. 

It will prove useful to introduce "reduced" enthalpy 
background densities 

la = 4TTGa 2 (p a + Pa ) , (All) 
7 = 4irGa 2 {p + p) = J2"/a ■ (A12) 

a 

From definition (A10) and Friedman equation (A4), 

7a=x a7 , 7 = 3(l + W )(W 2 +A')/2 . (A13) 

Finally, we give the rate of change of some of the above 
quantities with respect to the conformal time r. Energy 



conservation requires 

Pa = -3H{pa +Pa) , 

This and the Friedmann equations give 



H = 



1 d 2 a 



4irGa 2 



{p + 'ip) 



3 

— — (H 2 + K, 



(A14) 



(A15) 



By eqs. (A9,A14), 

w = m{l + w)(w-c 2 w ) . 
By eqs. (A7,A9) and (A14), 

7o = -W(l + 3c2) 7o! 7 = -H(l + 3c^) 7 . 
From the last two equations we also see that 



= -3W(c 2 - c 2 w )x a . 



(A16) 
(A17) 

(A18) 



2. Matter perturbations 

We parameterize T£ v perturbations by the particle 
number overdensity 20 S a = Sn a /n a , the peculiar velocity 
vector v a , and the anisotropic stress perturbation U l J: 

T„°0 = ~(Pa+Sp a ), Sp a = {p a +Pa)S a , 

T°i = ( P a+Pa)v la , (A19) 



S) (Pa +S Pa ) + ( Pa + Pa) , II* * = 



In scalar modes the 3- vectors vt a and 3-tensors 14^ are 
derivatives of scalar velocity potentials u a , 



Via = -VjU Q , 

and anisotropic stress potentials 7r a 
3 f_,_ 1 



n:, 



V t V J --S}V 2 )7r a 



(A20) 



(A21) 



With the normalization (A21), V;Vj-n^ = V 4 7r a . The 
potentials 7r a are related to some of the alternative vari- 
ables used to describe anisotropic stress as a = — V 2 7r for 
the variable of Ref. [57] and 14 = -(3/2)(l + p/p)V 2 ir for 
Ref. [60]. 

For the perturbations in the total energy-momentum 
tensor T M „, parameterized analogously to eqs. (A19), 

S = '^2x a S a , U^^XaUa, 7T = ^ X a TT a , (A22) 
a a a 

where x a were defined in eq. (A10). 



20 Many authors use 6 a for the energy overdensity Spa/ pa = (1 + 
w a )S a , in our notations. 
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3. Metric and Gauges 

Under a reparameterization of spacetime coordinates 
(gauge transformation) 

t = t + 6t(t,t), x? = x l + 6x i {T,r) , (A23) 

the perturbation variables transform as 

S a = S a + SHSt , u a = u a — St, it a = n a . (A24) 

Note that the following quantities are gauge invariant: 
S a + 3Hu a , or 5 a - Sb, or u a - u . 

The perturbed metric of space-time can be parameter- 
ized as 



ds 2 = a 2 (r) [(-1 - 2A)dr 2 - 2B t drdx 1 + 

+ (1 - 2H L )dr 2 - 2H ij dx i dx j ] , 



(A25) 



where H\ = 0. For scalar modes, the 3- vector Bi and 
the 3-tensor can be written as spatial derivatives of 
scalar functions: 

B l = V, b , H}= (v% - i <5]V 2 ) x ■ (A26) 

The metric perturbations transform under the gauge 
transformation (A23) as 



.4 


= A — 


ST - HSt , 


(A27) 


Bi 


= Bi- 


- V^t + S±i , 


(A28) 


H L 


= H L 


+ HSt + \\7tSx 1 , 


(A29) 




= HI- 


f I- (V^Sxa + VjSx 1 ) - - S i i V k Sx k . 
2 3 J 


(A30) 



For scalar perturbations with 

Sx.i = V; SX 

the potentials b and \ m e q- (A26) transform as 
b = b - St + SX , x = X + SX . 



(A31) 



In the conformal Newtonian (longitudinal) gauge the 
gauge conditions on scalar perturbations are &W = 
and x = 0- For brevity, we refer to this gauge as the 
"Newtonian". Defining $ = and * = H ( L N) we 

arrive at the metric of eq. (1). 

In the synchronous gauge one sets A( s ' = and 

(s) 

B\ ' = 0. The observers who are at rest in the syn- 
chronous gauge are free falling in the gravitational field 
and their locally measured proper time sets the coordi- 
nate time. By eqs. (A31), the gauge transformation from 
the synchronous to the Newtonian gauge is St 



-x 



(*) 



SX 



-X^ s - Hence, from eqs. (A27, A29), the Newtonian 



potentials are related to the scalar metric perturbations 
in the synchronous gauge as 

* = X (s) + Hx [s \ * = #i S) - Hx {s) - \ V 2 x (s) - (A32) 

For the reverse transformation from the Newtonian to 
the synchronous gauge, by eqs. (A27,A31), St(x) and 
SX(x) are any functions satisfying 



St + HSt = $ 



SX = St . 



(A33) 



The initial values i5r(ri n ,r) and (5A(ri n ,r) can be chosen 
arbitrarily, corresponding to the residual gauge freedom 
of the synchronous gauge. The metric perturbations in 
the synchronous gauge are obtained as 



H 



(s) _ 



= ^ + H5t + iv 2 <5A , x (s) = 5X . (A34) 



In the spatially flat gauge the scalar perturbations of 
the spatial part of the metric are absent: = and 

X^' = 0. By eqs. (A29, A31), this gauge has no residual 
gauge freedom. It is obtained from the Newtonian gauge 
with St = —ty/H, SX = 0. In terms of the Newtonian po- 
tentials, the scalar metric perturbations in the spatially 
flat gauge are 



A {f) = $ + * 



b (f) = 



H 



(A35) 



The comoving gauge is defined as the gauge in which 
the scalar components of the total Tf vanishes: vS c ^ = 
0. By eq (A24), the transformation from Newtonian to 
the comoving gauge is achieved with St = u 1 --^. For 
the second gauge condition, fixing Sx l in eq. (A23), it 
is convenient to choose x = 0- Then Sx l = 0. By 
eq. (A29), the spatial curvature potential in the comoving 
gauge is related to the Newtonian gauge variables as 



Tl EE ttf = * + HU^ 



(A36) 



For the remaining scalar metric perturbations 
eqs. (A27, A31) give 



A (c) = § _ _ Hu (N) 



u^> . (A37) 



By the following eq. (A43) and eq. (A13), the "comoving 
curvature" 1Z can be easily transformed into its conven- 
tional form 




(A38) 



The uniform density gauge corresponds to the con- 
dition 5^ ee 0. Hence, eq (A24), it is obtained from 
the Newtonian gauge with St = —S^ N y(3H). Taking 
-^(") = o to be the second gauge condition, one finds 
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that the 3-curvaturc potential Hr, in the uniform density 
gauge is 



C = = * - I 8™ 



(A39) 



By Einstein equations (A42-A43), the curvature pertur- 
bations in the comoving, uniform density, and Newtonian 
gauges are related as 



V 2 * 



(A40) 



4. Conservation and Einstein Equations 

Newtonian gauge: 

The dynamics of the species density and velocity per- 
turbations follows from the conservation law T^ v - v = 
as 



momentum conservation equations give 



5 { a ] = V 2 u { a s) + 3H { L S) , 



= c 2^)_ XaU (s) +V 2^ , 



(A48) 



The corresponding linearized Einstein equations, e.g. 
Ref. [57], in spatially flat background are 

V 2 (h[ s) - ^V 2 x (s) ) - 3HH ( L S) = 7 <S (S) , (A49) 

H ( L S) - iv 2 x (s) = 7" (s) , (A50) 
ffW +W iiW= 7 (| + i)^), (A51) 

X (s) + 2Wx (s) - (h^ ] \^X { A = - 3 7 7T . (A52) 



5. Dynamics in phase space 



5a = V 2 U Q + 3* , 

Ua = C 2 (5 a - XaU a + V 2 7T a + $ 



(A41) 



(for the scalar mode), where Xa = W (l — 3c 2 ) is the 
Hubble drag rate for the species a. The evolution of 
the anisotropic stress potential ir a is determined by the 
internal dynamics of the species. 

The linearized Einstein equations in Newtonian gauge, 
Ref. [57, 99], are easily reduced to 



(A42) 
(A43) 



V 2 #-3W(* + W$) = jS , 
* + HQ = 7« , 



Sp 



* + n (2* + $) - 'swh 2 ® = i-f p 5 + 1^ 2 k, 



* - $ = 377T , 



(A45) 



where 7 is introduced by eqs. (A12,A13) and the back- 
ground is assumed spatially flat. If for all the species 
Sp a = c-l^pa then (Sp/Sp)S = J2a x a, c l^a- By equa- 
tions (A42-A43), 



V 2 * = 7<S (c) , 



(A46) 



where 8^ = ^2 a x a 8 a c ^ is the averaged particle number 
density perturbation in the comoving gauge, 



5^ = s a + m% 



(A47) 



Six variables specify the coordinates of a particle in 
phase space at a given time. For them, following Refs. [3, 
57], we take the comoving coordinates of the particle r l 
and the comoving momenta 



apt 



(A53) 



where pi are the proper momenta measured by a comov- 
ing observer, who is at rest with respect to the coor- 
dinate frame. For a particle with a mass m a the mo- 
menta that are canonically conjugate to the variables r % 
are Pi = m a dxi/ V 1 —ds 2 = (1 — i/Oft- 

The particle density in phase space is specified by the 
canonical phase space distribution f a (r l , Pj , t): 



dN a = f a (r\P j ,T)d d r- l d d P j 



(A54) 



for every species of particles and their states of polariza- 
tion a. The energy-momentum tensor of the species a is 
given in the Newtonian gauge by the following simple ex- 
pression up to the first order of cosmological perturbation 
theory [3, 57]: 



(A55) 



with p° ee — Po ee y/ (q/a) 2 + ml and p l = Pi = q x /a. 
Below we drop the species label a when referring to any 
sort of particles in general. 

The evolution of the phase space distributions obeys 
the Boltzmann equation: 



Synchronous gauge: 



dr l dq ' ' drii \dr 



(A56) 



In Appendix B we refer to the evolution equations in 
the synchronous gauge. In this gauge the energy and 



where / is considered as a function of the coordinates 
r l , q ee |<fo|, n% ee qi/q, and r. The right hand side of 
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eq. (A56) describes the change of the phase space den- 
sity due to particle collisions. Linearizing the Boltzmann 
equation relative to an unperturbed background phase 
space distribution, / = / + Sf, in the Newtonian gauge 
one finds, Ref. [57]: 



(Sf)' + 7 V, (Sf) + q J£ (it - - q niVi$ 



(A57) 



where e = \J q 2 + a 2 m 2 . 



APPENDIX B: LOCALITY OF ADIABATIC 
GREEN'S FUNCTIONS 

In this Appendix for any "reasonable" cosmological 
system we construct the perturbation Green's functions 
in that all the scalar gravitational and matter distribu- 
tion potentials are identically zero in the Newtonian or 
synchronous gauges beyond the Green's function particle 
horizon \x\ > t. We argue that the perturbations formed 
with these Green's functions by their convolution with 
any smooth kernel are adiabatic, as defined in footnote 8 
on p. 6. A Green's function Fourier component with any 
wave number k presents a non-zero non-decaying adia- 
batic mode. 

We assume that all the scalar perturbations of mat- 
ter distributions, classical fields, or the metric tensor in 
our system can be parameterized by a sum of 3-scalar 
functions (matter potentials) acted by polynomials of V*. 
This is easily achieved for classical fields that are scalars, 
vectors, or tensors of a higher integer rank. Complica- 
tions would arise for spinor and other non-integer spin 
fields. However, being fcrmions, such fields can not sup- 
port coherent semiclassical excitations. Without sacri- 
fice of generality, we require that the polynomials of V, 
are homogeneous I = 0, 1, etc. degree polynomials that 
transform under spatial rotations as irreducible tensors 
of the rank I. We treat all the dynamical degrees of 
freedom, whether they are dynamical "coordinates" or 
"velocities" , on equal footing, implying the Hamiltonian 
formulation of the system classical dynamics in the coor- 
dinate time of the chosen gauge. 

We set the Green's function initial conditions at a spa- 
tial slice of an infinitesimal thickness at a coordinate time 
ri n > 0, which is eventually sent to zero. Within this 
slice we impose the Newtonian gauge conditions 6=0 
and x = on the scalar metric potentials (A26). We 
require that, in the Newtonian gauge, all the matter dis- 
tribution or field perturbation potentials initially vanish 
for |a;| > r; n . (The Green's functions considered in this 
paper are homogeneous in the normal to x spatial direc- 
tions y and z, Sec. IV.) Inside the interval \x\ < T m , 
to be specific, we set all the I > 1 matter distribution 
potentials at T m to zero and adjust the initial conditions 
for classical fields, if any are present, so that u(Tj n ) = 0. 



(Then the metric on the initial spatial slice satisfies the 
comoving gauge condition as well. The spatial slices in 
the Newtonian and comoving gauges will, in general, dif- 
fer for t > Tj n .) At the time T m by eq. (23) 



V 2 * - 37* = 



(Bl) 



Let us demand that 9(r m ,x) = (l/r; n ) ^>o(x/Ti n ) 
where 'I'o(x) is an arbitrarily chosen even function that 
vanishes for \\\ > 1 and that has a non-zero integral 
over dx from —1 to 1. Then the coordinate "particle 
number" density perturbation c?(rj n , x) should be given 
by the right hand side of eq. (Bl). The corresponding 
proper number density perturbation S at r; n is 



S = d + 3* 



V 2 * . 



7 



(B2) 



We set all the species initial densities and the initial con- 
ditions for classical fields so that S a (ri n ,x) = (1/7)V 2 ^ / 
and all the other matter potentials arc even and vanish 
for | a; | > Tin- For interacting species, the separation of 
the total energy-momentum tensor T^ v into T£ v and so 
the definition of the species density perturbations S a may 
be ambiguous. In this case, we take any of the possible 
definitions. In the limit T m — ► all of them lead to the 
physically identical adiabatic Green's functions. If any 
of the matter potentials or fields remain unfixed, we take 
them unperturbed at Tj n for all x. 

The motivation for these initial conditions is to con- 
struct an initially localized almost "pure" curvature per- 
turbation. For any Tj n > 0, the curvature perturbation 
must be generated by some non-zero matter disturbance. 
The above procedure attempts to restrict the required 
matter inhomogeneity so that in the limit T- m /\ — > 
only curvature appears to be perturbed. Indeed, after 
the convolution of the constructed Green's function with 
a primordial fluctuation field A(x) that is smooth on the 
comoving scales below some A > 0, the number density 
perturbation of any species at r; n by eq. (Bl) tends to 



lim (S(x)) = lim ( dx'A(x') Y-*iZ — ^1 
ri n — >o n n —>o J 7 

1 dxMx)) lm:^ 

j J r in -0 7 



= (B3) 



The last limit vanishes as T 2 n /X 2 . Since the other ini- 
tial matter distributions are chosen unperturbed, they 
remain unperturbed after the smoothing. However, the 
Tm — * limits of the smoothed comoving curvature per- 
turbation 1Z = * + Ttu as well as the variable £ = — d/3 
at r; n are non-zero: 



lim o (K(x))=lim o (ax))=U d x *o(x)JA(s) . (B4) 

Given the locality of the Green's function, which is 
proved next, and the initially vanishing smoothed matter 
perturbations, we argue at the end of Sec. II B that the 
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smoothed perturbations £ and 1Z remain constant while 
r«A. The Fourier transformation is a convolution with 
the kernel A(x) — exp(— ikx). Therefore, it gives non- 
decaying non-zero curvature perturbation modes (B4), 
as long as J_ 1 d\ is chosen different from zero. 

The smoothed density perturbation (S(x)) will remain 
as small as 0(t 2 /A 2 ) for r- m < t < A in the comoving 
gauge, where it is always related to <3> by eq. (A46). In 
other gauges, including the Newtonian, the proper den- 
sity perturbation (S(x)) and other proper matter pertur- 
bations do not need to vanish for T m < r <C A. For 
example, the species proper density and velocity pertur- 
bations are given by eq. (A24), where St is the time lapse 
between the comoving and the considered gauge. 

Now we prove that in a synchronous gauge and the set 
initial conditions at T ln all the matter and gravitational 
Green's functions remain zero for r > T ln beyond the par- 
ticle horizon |x| > r. The coordinate transformations St 
and Sxi = V,<5A from the Newtonian to the synchronous 
gauge can be chosen so that <5r(T in ) = and 5\{T m ) = 
for all x. Then, from eqs. (A33), Sf(Ti n ) = $(ri n ) and 
^A(ri n ) = 0. Eq. (A34) gives that the synchronous gauge 

potentials and x^ and their first time derivatives 
vanish at Ti n for \x\ > T m . Same is true for the trans- 
formed initial perturbations of all the matter distribu- 
tions and their rates of change. 

Causality requires that in the synchronous gauge the 
matter and the metric remain unperturbed beyond the 
particle horizon of the original perturbation. This means 

that H ( L S \ V 2 x (s \ c.f. eq. (A26), Vu( s >, V 2 tt, etc. are 
zero for |x| > r. We argue that the potentials \ \ u ^ s \ 
7r, etc. themselves are also zero for \x\ > r. As an illustra- 
tion, let us consider the Euler's equation for v x = —V x ii 
that in the synchronous gauge reads 



Vr u 



V, \c 2 J^ - H(l - 3c2>« + V 2 tt1 . (B5) 



The gradient in this and the similar equations can be 
dropped after we define the expressions under the gradi- 
ents to be equal at a certain xq, which may depend on r. 
Consistent definition of the potentials is obtained if Xq 
is chosen outside of the perturbation horizon; suppose, 
xq(t) < — t. Since a sufficiently high spatial derivative 
of every synchronous gauge potential gives a perturbation 
of some matter distribution or a metric tensor component 
and thus, by above, vanishes for x < — t, the potentials 
themselves can be chosen to vanish in that interval for 
all t. We prove that the potentials then also vanish for 
x > T. 

Let {/;} be multipolc potentials for some matter distri- 
bution perturbation F(r, n), which may depend on other, 
not displayed, phase space or internal coordinates, 



F(r,n)^£(-l)'(2Z + l)P i 



riiV, 



V'/i , (B6) 



where i runs over the spatial coordinates x, y, and z. We 
suppose that the distribution perturbation dynamics can 



be described as 
F(r,n) = 



/ local in space functional 
V of perturbation fields 



(B7) 



The functional on the right hand side is linear in the lin- 
ear order of perturbation theory. Multiplying both sides 
of eq. (B7) by a spherical function Yi m (n), integrating 
over the solid angle d 2 f^, and applying the identity 



J d 2 Q A Y lm (n) Pi 



27TIM?) 



on the left hand side of eq. (B7) we find /; times a ho- 
mogeneous n-degrec polynomial Qi m (V x , V y , V z ). The 
polynomials {Qim}' n= _/ transform under spatial rota- 
tions as Yi m . Since all the perturbation potentials are 
defined to be 3-scalars, Vj are the only quantities that 
may appear on the right hand side and that transform 
non-trivially under spatial rotations. Therefore, after the 
convolution with Y; m (n), every term on the right hand 
side of the linear equation (B7) must contain the same 
polynomial Qi m (^x, Vy, V z ) as that appears on the left 
hand side. Putting m to and applying the resulting 
evolution equation to the y and z independent Green's 
functions, we find (V x ) 1 fi on the left and at least I deriva- 
tives Vx in every term on the right hand side. After I 
integrations over dx with zero initial values at x < — r 
we find an equation of the form 



j /local in space functional \ 

\ of scalar potentials J 



(B9) 



The corresponding equations for the studied in the 
main text system of photon-baryon and CDM fluids 
and ultra-relativistic neutrinos are eqs. (4-5,20). The 
scalar potential evolution can be generalized to describe 
photon-baryon Thompson scattering and photon polar- 
ization as shown in the Appendix of Ref. [56] . The evo- 
lution of the synchronous metric potentials \ x^ s \ 
and can as weu b° presented in the form (B9) using 
eqs. (A49,A50, A52). 

Thus the dynamics of all the potentials in the syn- 
chronous gauge can be reduced to an initial value 
Cauchy's problem. The initial conditions at r = Tj n 
were chosen even. Assuming the dynamical equations 
are invariant under x — > —x inversion, the resulting solu- 
tion of the Cauchy's problem will remain even for all r. 
Hence, the potential values for x > r are equal to those 
at x < — r, i.e., vanish in the synchronous gauge. 

The reverse transformation to the Newtonian gauge is 
achieved with St = ~x and SX = — x j see Sec. A3. 
These functions have been proved to vanish outside of the 
horizon. Therefore, all the matter multipole potentials in 
the Newtonian gauge and the gravitational potentials <I> 
and ty, related to the synchronous metric perturbations 
by eqs. (A32), as well vanish for \x\ > r. 
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APPENDIX C: 0(R„) ORDER CALCULATIONS IN THE RADIATION ERA 

Gravitational potentials 

In Sees. IVD-IVE we derived integro-diffcrcntial equations satisfied by the potentials $± = ± $)/2 during 
radiation domination. We can rewrite these equations, eqs. (101 ,100) and (106, 105, 114-115) as 



2CinX(l-X 2 )^(l-|xl)-6 / d X ' I d X 



n $ + (x")x'P2(x') + 3>+(x')x"P2(x") 



x" - x' 



and 



with 



Mx) 



F-(X) 



P<s> 



V3 



F-(X) 



dx' 



,12 _ I % 



■*'_(x') 



V3(l - 
1-2FL 



dx^-(x) 



(Cl) 

(C2) 

(C3) 
(C4) 



In this Appendix we use them to calculate in 0(R V ) order the prefactors at the singular terms in the photon and 
CDM density Green's functions. Up to small 0(R%) corrections, this yields the density modes on subhorizon scales. 

The quantity in 0(R U ) order is obtained by substituting the zeroth order solution for $ + , eq. (109), into the 
right hand side of eq. (Cl). The substitution gives 



<!>'_ 



iUin { Mx) - h signx] 0(1 - \ X \) - Mx) - h signx] 



V3 



0(R 2 U ) 



where 



(3x 2 - I) 3 



In 



V3 



Hx) 



3\/3 



(i - x 2 ) 



4V3 

2x (3 X 2 + 1) -(3 X 4 + l)ln 



X 7% 



i + x 



(C5) 



i-x 



Jo = Hi) = J 2 (^= 



2 - — ln(2 + V3) 



Photon singularities 

The magnitude of the photon density acoustic spike p 7 , 
eq. (115), can be calculated analytically using F- defini- 
tion (C3) and the above expression for <&'_. For this we 
employed the symbolic calculation program "Mathcmat- 
ica" (http://www.wolfram.com). The resulting expres- 
sion turns out rather lengthy but is easily evaluated to 

Vl ~ | Cm [1 - 0.2683^ + 0{Rl)\ . (C6) 

The residue r 7 , as defined by eq. (118), follows from 
eqs. (112,106) as 

ri = $ + (-L) =( in R v (A _ i\ -0.1656 ( m R v . (C7) 



Substitution of the found values in eq. (120) leads to the 
results (121-122). 



CDM singularities 

Now we proceed to the calculation of the neutrino 
corrections A c and 8c in the CDM density perturba- 
tion (128). The equation for the radiation era CDM 
Green's function (125) gives 

(Xdc)' = ~$'+PMX) ■ (C8) 
X 

Integrating this expression, fixing p c \ by the requirement 
that xdc is an odd function that vanishes at |x| > 1, and 
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dividing the result by x, we obtain 
G( X ) G(l) 



d c ■ 
where 



X 



G(x) 



0[l-\x\)+PcSu{x) , 



dx $'(x) 

X 



(C9) 



(CIO) 



The initial condition J d\ d c = — 3Cin fixes the new 
constant p c as 



Pc 



-3Cin - / d\ 



G(x) G(l) 



X 



The dx integral of G(l)/|x|, which would be infinite in 
Riemann sense, is equal to in the sense of generalized 
function integration, see Table I on p. 15. Hence, 



Pc 



-3Ch 



1 d X G(x) 



X 



(Cll) 



The CDM density perturbation modes in the radiation 
era are given by the Fourier components of the Green's 
function (C9). In the radiation era for adiabatic initial 
conditions the potential 3>(x), generated by the photon 
and neutrino perturbations, is regular at x = and even. 
Hence, the function G{x) in eq. (CIO) is regular at X = 
and odd, and so is G(x)/x regular at the origin. The 
subhorizon (kr 3> 1) values of the CDM Fourier modes 
are fully specified by the singular terms in eq. (C9), which 
are proportional to l/\x\ and <5d(x)- From Table II on 
p. 16 we find that 



d c (r, k) = A c (In <p s + c) + 0(^ _1 ) 



with ip s = kr/y/o and 



In 3 



A C = 2G(1), c = ~i+ — + 



Pc 



(C13) 



2 ' 2G(1) 

In the R v -> limit, with the potentials (108-109), 



-3Ci, 



X 



V3 , 



signx , \x\> 4 



Then = -3C in and 



= -3Cin 



(C15) 
PcSuix) , 



where by eq. (Cll) 



=3Gn(l + ln3) . 



(C16) 



Fourier transforming eq. (C15-C16), where the singular 
function l/|x| can be transformed using Table II, we 
obtain the CDM density perturbation mode (127). Of 
course, the radiation era, R v — * Fourier mode could be 
obtained directly in fc-space by integrating the evolution 
equation (124). Now we consider how the CDM density 
perturbation changes when neutrinos are added. 

Although the analytical calculation of the integrals 
in eqs. (CIO, Cll) in the 0(R V ) order may be possi- 
ble, it does not appear easy. On the other hand the 
numerical evaluation of the absolutely convergent inte- 
grals in eqs. (CIO, Cll), given the potentials (C5, 106), is 
straightforward and yields 



GO) 

Pc 



-3Cin [I + 0.2297 R u + 0(Rl)] , 
3Ci„ [1 + In 3 + 1.746i?„ + 0{Rl)] 



(C17) 
(C18) 



(C12) Then for A c and c, eqs. (C12-C13), 



A c ~ -6Cin [1 + 0.2297^ + 0^) 
c ~ 7 - | - 0.6323^ + 0(Rl) . 



The corresponding values for A c = A c /A 



(fl„-o) 



(C19) 



1 and 



(C14) 6c = c 



are used in Sec. IV F, eq. (129). 
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